Skillful prediction of tropical Pacific fisheries provided by Atlantic Niños

Annual catch data from the Food and Agriculture Organization of the United Nations (FAO) major fishing areas Eastern Central Pacific (#77) and Southeastern Pacific (#87) for the period 1950–2014 were retrieved from the Sea Around Us Project (SAUP) database (www.seaaroundus.org/data/#/fao). These are a combination of official reported data from FAO FishStat database (www.fao.org/fishery/statistics/en) and reconstructed estimates of unreported data [43]. All types of marine species and fishing methods were considered altogether in the SAUP timeseries, including catches from exclusive economic zones. To retain interannual variability and minimize the potential impact of external factors (e.g. fishing effort, climate change impact) on historical catch data, a high-pass Butterworth filter with 11 year cut-off period was applied to the time series.

Observational 1° longitude-latitude SST data were retrieved from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST [44]). Simulated physical and biogeochemical data for the historical period (1971–2004) were obtained from the Geophysical Fluid Dynamics Laboratory (GFDL) Carbon Ocean Biogeochemistry and Lower Trophics run (COBALT [45]). Physical-biological parameters from this simulation (so-called GFDL ocean reanalysis within the FishMIP community) and additional datasets are fully described in table 1. The COBALT ocean-ice simulation, performed with the Modular Ocean Model v4p1 [46], has a 1° horizontal resolution (1/3° along the equator) and 50 vertical layers with 10 m increase in the top 200 m. The simulation is forced by air-sea fluxes from the CORE-II data set [47], a blend of the National Centers for Environmental Prediction (NCEP) atmospheric reanalysis [48] and satellite data. Whereas sea surface salinity is restored to observations in GFDL-COBALT due to drift issues, SST and biogeochemical outputs are free running (hindcasts [45, 49, 50]). Full information of the GFDL-COBALT validation against physical-biological observations (e.g. mixed layer depth, nutrient concentration, phytoplankton/zooplankton biomass, carbon/energy fluxes) is available in Stock et al [45].

In line with CORE-II forcing on GFDL-COBALT, NCEP reanalysis was selected to analyze atmospheric circulation dynamics and variability in this work (cf table 1).

Model simulations from the FishMIP v1 intercomparison protocol [42] were analyzed. Available model outputs through the ISIMIP repository (https://esg.pik-potsdam.de/search/isimip/) were considered for the analysis: EcoOcean, BOATS and Macroecological [51].

EcoOcean is a complex spatial-temporal global ecosystem model that comprehensively simulates interactions within all trophic levels and species groups using the Ecospace [52] spatial-temporal model at its core. The depth dimension is implicitly considered in the model by preferent habitat patterns and food-web interactions (one dimensional). EcoOcean v1 uses SAUP effort data [53] spatially distributed across large marine ecosystem regions as fishing forcing. Effort data combine fleet information (e.g. number and type of boats) and dedicated time for fishing. The forcing is applied as a mortality rate for functional groups in the simulation through the removal of existing biomass per unit of time. The model is able to replicate historical catch successfully on a global scale. A full description of EcoOcean is available in Christensen et al [54].

BOATS is a global ecological model that determines size spectra of fish biomass as a function of net primary production and local temperature. Biomass production in higher trophic levels is limited by available photosynthetic energy, temperature-dependent growth rates and trophic efficiency scaling. The model enables a realistic representation of biological and ecological processes, allowing evaluation of the links between climate, biogeochemistry and upper trophic dynamics. BOATS employs a bioeconomic approach based on SAUP catch price data to determine spatial-temporal changes in fishing efforts. The vertical dimension is a single surface-integrated layer and movement of fish is disregarded. Full information on the model architecture is available in Carozza et al [55].

Macroecological is a static-equilibrium model that relies on predator-prey mass ratios, transfer efficiency and metabolic demands dependent on body mass and temperature to predict size and distribution of marine consumers (as a whole) at any given time and location. It has a single vertical (surface-integrated) layer and movement of fish is disregarded. Extended information on the model is available in Jennings and Collingridge [56].

Seven different simulations were analyzed in this study. These include fished and unfished runs with and without diazotrophs dynamics (table 2). All historical runs (1971–2004) are forced by GFDL-COBALT simulation, whereas an EcoOcean 2021–2054 projection is driven by IPSL-CM5A-LR RCP8.5 [57]. In the latter, fishing forcing is kept constant at 2005 rates [42]. Although GFDL-ESM2M was also available as driving forcing for 2021–2054, it was eventually disregarded due to its inability to simulate a connection reminiscent of the observed Atlantic-Pacific teleconnection subject of this study [58]. Complete information on simulation settings, acronyms and data sources are available in table 2. All common FishMIP model output variables considered are detailed in table 1 and are publicly available through the FishMIP global model repository [51]. The most detailed model simulation (EcoOcean) forced by historical fisheries effort and which includes diazotrophs dynamics is used as reference throughout the manuscript (hereafter EcoG-Fis), as it provides the closest possible approximation to reality.

Monthly anomalies of physical-biological variables were calculated by removing their total period average to the time series of each grid point. Seasonal (4 month) anomalies were computed as aggregated monthly anomalies, with linear trends effectively removed to minimize the impact of anthropogenic climate change [59].

Physical-biological indices were determined by spatially averaging seasonal anomalies over regions of interest, e.g. Niño3 (150° W–90° W, 5° S–5° N) and Atl3 (20° W–0° E, 3° S–3° N). All indices were normalized by dividing the time series by their standard deviation.

Whereas Atl3 was computed for SST on boreal summer (July–September; JJAS) to represent Atlantic Niño variability during its peak intensity [22] (a.k.a. Atlantic Equatorial Mode), Niño3 and global seasonal anomalies of physical-biological variables (cf table 1) were calculated as rolling 4 month averages between monthly lags −8 to 24, where lag 0 corresponds to JJAS, lag 1 (−1) to JASO (MJJA) of the same calendar year and so forth (back).

Statistical significance for correlation outputs (figures 1, 3 and 4; 95% confidence interval) was evaluated using a Student-t test that accounts for the autocorrelation of the time series and the calculation of effective degrees of freedom [60].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The SST-based Statistical Seasonal forecast model (S4CAST) is a novel statistical seasonal model that utilizes SST data (predictor) to forecast the response of any given impact variable (predictand). To this end, the model performs a maximum covariance analysis (MCA [61]) between a predictor Y field (SST anomalies) and a predictand Z field for any given region, time period and lag.

MCA is based on single value decomposition applied to the non-square covariance matrix (C) of both Y and Z fields. The method calculates linear combinations of the time series of Y and Z (a.k.a. expansion coefficients U and V, respectively) to maximize C.

A leave-one-out cross-validated hindcast based on MCA outputs is then performed [62]. For this purpose, predictor/predictand data combinations from all years except the year being predicted at each training stage are selected for model training and validation [62, 63]. A Monte Carlo test with 1000 random permutations is applied to MCA (U vs. anomalous field correlations) and cross-validation of skill scores for hypothesis testing (95% confidence interval). The Pearson correlation coefficient is selected as a qualitative indicator for model skill (sign of response).

Further details on the S4CAST statistical model framework are available in Suárez-Moreno and Rodríguez-Fonseca [64].

We find year-to-year variations of historical global SSTs and eastern Pacific total catch records from the SAUP database (FAO areas 77 and 87) robustly correlated for the period 1950–2014. The years when observed fish catches are higher than usual synchronize with well-developed La Niña-like SST anomalies over the tropical Pacific and elsewhere [11, 65] (figure 1(A)). Lagged SST anomalies from the previous calendar year depict a preceding Atlantic Niño signal, together with a developing La Niña over the Pacific (figure 1(B)).

Likewise, we consider EcoOcean model historical simulation in a fishing forcing scenario (EcoG-Fis) for the same analysis. We first characterize the tropical Pacific in terms of biomass density of commercial species. As expected, this is the most productive and fluctuating region worldwide in the model (figures 1(C) and (D)). Despite the slightly smaller area and shorter period considered (FishMIP coverage; 1971–2004), results very similar to observations are obtained from the model (figures 1(C) and (D)). Apart from equatorial Atlantic SST anomalies and intrinsic 2–7 years Niña-Niño intermittency, no other significant relations with tropical SST anomalies (e.g. NTA, Indian) are found, either in observations or in simulations during the years preceding annual-averaged total catch anomalies (figure 1 and supplementary material figure S1 (available online at stacks.iop.org/ERL/16/054066/mmedia)).

Correlation between the Atlantic Niño index (Atl3 JJAS SST; hereafter Atl3-Had) and the tropical Pacific total catch corroborates the lagged one year Atlantic-Pacific relationship in the model and observations (supplementary material figure S2). As expected, this relationship substantially weakens and changes sign in subsequent year lags. In agreement with Rodríguez-Fonseca et al [16], a slightly stronger Atlantic-Pacific connection is found for 1971–2004 compared to 1950–2014.

Owing to constraints in historical fishing records (averaged annual data over vast areas), we focus on gridded 1-degree monthly FishMIP modeling outputs for fisheries prediction. Again, we consider the most detailed model forced by historical fisheries efforts (EcoG-Fis), which provides the closest approximation to reality (table 2). We then use the Sea Surface temperature-based Statistical Seasonal forecast model (S4CAST [64]) to perform a cross-validated hindcast of simulated tropical Pacific marine yields: total catch at sea (tc) and total system carbon biomass density (tsb).

Consistent with figure 1, we use as predictor (Y) equatorial Atlantic SST anomalies during the peak intensity season (Atl3 JJAS-lag 0) of Atlantic Niños/Niñas [16, 22, 66] on GFDL-COBALT, the driving historical forcing of FishMIP. Tropical Pacific (150° W–75° W, 20° S–20° N) simulated total catch and total system biomass anomalies are used as predictands (Z) for different monthly lagged periods (0 to 42).

For all monthly lags, Atlantic Niño emerges as the main predictor of tropical Pacific total catch anomalies, which can be hindcasted over large-scale areas and up to three years in advance (figures 2(A)–(E) and supplementary material figure S3). The quantification of the surface area of the domain with skillful prediction returns values of up to 3.4 million km² at lag 20 (figure 2(F)). If total system carbon biomass density is considered as predictand (Z) instead of total catch, final results are very similar but shifted ∼9–12 months back in time (compare supplementary material figures S3 and S4), with a maximum of 4.2 million km² of surface area with skillful prediction at lag 8 (figure 2(F)).

To establish a physical-biological mechanism explaining the results above, we evaluate tropical Atlantic (Atl3) SST impact on Pacific (Niño3) region for different area-averaged physical-biological variables and monthly lags. The simulation selected for this purpose is again EcoG-Fis, although other marine ecosystem models (BOATS, Macroecological) and scenarios are also considered (cf tables 1, 2 and Tittensor et al [42]).

First, we check whether the Atlantic-Pacific physical teleconnection is realistically represented in GFDL-COBALT. This is motivated by the fact that SSTs are not restored to observations in the simulation and equatorial wind fields in CORE-II depict inaccuracies [50, 67]. Particularly, inconsistencies in wind speed, humidity, and solar radiation have been reported between NCEP and CORE-II [47]. In spite of this, the good performance of GFDL-COBALT representing Atlantic Niño variability and associated impacts is confirmed in figure 3(A) and supplementary material figure S5, where results obtained are very similar to observational studies [16, 28]. For a warming in the equatorial Atlantic, a strengthening of the Pacific Walker Circulation is observed, together with increased subsidence and anomalous easterlies over the western equatorial Pacific, features that promote subsequent La Niña conditions (supplementary material figure S5). Similarly, previous studies have shown that anomalous easterly winds over the western equatorial Pacific promote the propagation of equatorial upwelling Kelvin waves [23, 31, 32, 68], thus favoring the emergence of nutrients from the seabed.

The linear correlation value between Atlantic Niño indices (Atl3 JJAS 1971–2004) in GFDL-COBALT and HadISST is 0.87, whereas that for Niño3 (NDJF) is 0.96 (p-values < 0.05). The fact that both SST variabilities are significantly governed by the movement of water in response to wind variations may be behind the good agreement between GFDL-COBALT and HadISST [39, 49].

Lagged-correlations for biological variables in figure 3(A) enable characterization of the bottom-up propagation of Atl3-related signal (hereafter Atl3-GF) across the Niño3 region food web. Closely synchronous to the growth of cold nutrient-rich SST anomalies (figure 3(A)), large phytoplankton production (lphy), total system (tsb) and total consumer (tcb) carbon biomass densities grow rapidly before monthly lag 0, and decay thereafter at different rates. The different declining rates are due to distinct variable definitions and growth rates of species included in each category. While lphy is large-phytoplankton dependent, tcb is dominated by zooplankton and tsb is a blend of both (cf table 1). Consistent with previous studies [69], we find no SST response by small phytoplankton (not shown). Variables negative lag behavior can be partially explained by construction of Atl3 index at its peak season instead of its growing stage (supplementary material figure S6), and the fact that ocean chlorophyll tends to slightly precede ENSO-related SST anomalies due to wind-iron supply dynamics [70].

Results for carbon biomass density of fish consumers greater than 10 cm (b10) and 30 cm (b30) depict a bottom-up propagation of the physical forcing over time, peaking ∼9–12 months later than large phytoplankton productivity (figure 3(A)). By construction of b10/b30, where fish species are classified in terms of body size, bottom-up indirect effects of climate anomalies seem to overbalance species-specific direct temperature-recruitment processes [65]. A feasible process for this bottom-up propagation is a recruitment increase of short-lived species caused by enhanced food supply at lower trophic levels [36].

The spatial response of physical-biological variables for lags 0, 9 and 24 is shown in figures 3(B)–(D). At lag 0, cold SST anomalies and enhanced lphy/tsb are present over the Niño3 region (figure 3(B)). At lag 9, a mature La Niña is present, together with increased b10 over the Niño3 area and lphy/tsb to the north and south. This northward-southward advection of lphy appears related to Ekman transport induced by enhanced trade winds (supplementary material figure S7). At lag 24 increased b10 also appears displaced to the north and south, while b30 is only located to the north.

A broader picture of this evolution is given in figure 4 in the form of Hovmöller diagrams. Enhanced lphy/tsb/tcb tend to overlap with cold SST anomalies, preserving the timing from figure 3(A). Growth starts at small negative lags over the central tropical Pacific, then anomalies propagate to the east until lag 3 (figure 4(A)) and split thereafter into two branches to the north and south (figure 4(B)). The long persistence of lphy/tsb/tcb anomalies may be caused by enduring presence of nutrients (nitrate and dissolved iron). As opposed to SSTs, surface nutrients do not interact with the atmosphere [8, 36]. Lastly, whereas b10 show a similar response to lphy ∼9 months later, increased b30 are constrained over the northeast tropical Pacific from lags 14 to 24 (figures 3(D) and 4(B)). The latter may be caused by regional currents disposition leading to latitudinal asymmetries in large fish abundance (cf contours in figure 1(C) as a hint). Nevertheless, a species-specific analysis (beyond the scope of this letter) is needed to answer this question.

The evolution of physical-biological variables shown in figures 3 and 4 are in line with the shape and progression of total catch and total system biomass predictable regions in figure 2 (large crosses) and supplementary material figures S3 and S4. This is consistent with the propagation in time of climate anomalies within the food web.

To evaluate the role of diazotrophs dynamics in the physical-biological mechanism, the analysis in figure 4 is replicated in the EcoG-fis-wodiaz simulation (identical to EcoG-fis but without diazotrophs dynamics, cf table 2). The differences obtained between the two simulations are exiguous and therefore not shown. This seems consistent with the fact that primary production linked to nitrogen fixation is much less compared to total primary production, particularly in highly-productive marine areas [42, 71, 72].

To test the sensitivity of ecosystem response to top-down fishing control and model setup, in supplementary material figures S8–S10 we perform the same analyses on EcoOcean and BOATS fished and unfished model output scenarios (table 2). Results indicate that climate-driven bottom-up propagation dominates the response throughout the trophic web, while fishing visibly affects large predators over large-scale areas. The differences observed between models (less persistent and westward displaced anomalies of all-size fish biomasses in BOATS model outputs) may be caused by different model complexity, internal hypothesis regarding species dynamics and species relationships, ecosystem functioning and fisheries schemes [42, 59].

Lastly, a similar analysis is conducted for Macroecological, although it only provides annual model outputs, does not include a fishing scenario, and propagation in time of anomalies through the food chain is unresolved (supplementary material figure S11). As expected, the model returns the same picture as in figures 3(B) and (C) for lagged years 0 and 1, with the difference that consumer response (tcb/b10/b30) to the forcing is no longer lagged in time.