Extreme and compound ocean events are key drivers of projected low pelagic fish biomass

Ocean extreme events, such as marine heatwaves, can have harmful impacts on marine ecosystems. Understanding the risks posed by such extreme events is key to develop strategies to predict and mitigate their effects. However, the underlying ocean conditions driving severe impacts on marine ecosystems are complex and often unknown as risks to marine ecosystems arise not only from hazards but also from the interactions between hazards, exposure and vulnerability. Marine ecosystems may not be impacted by extreme events in single drivers but rather by the compounding effects of moderate ocean anomalies. Here, we employ an ensemble climate‐impact modeling approach that combines a global marine fish model with output from a large ensemble simulation of an Earth system model, to identify the key ocean ecosystem drivers associated with the most severe impacts on the total biomass of 326 pelagic fish species. We show that low net primary productivity is the most influential driver of extremely low fish biomass over 68% of the ocean area considered by the model, especially in the subtropics and the mid‐latitudes, followed by high temperature and low oxygen in the eastern equatorial Pacific and the high latitudes. Severe biomass loss is generally driven by extreme anomalies in at least one ocean ecosystem driver, except in the tropics, where a combination of moderate ocean anomalies is sufficient to drive extreme impacts. Single moderate anomalies never drive extremely low fish biomass. Compound events with either moderate or extreme ocean conditions are a necessary condition for extremely low fish biomass over 78% of the global ocean, and compound events with at least one extreme variable are a necessary condition over 61% of the global ocean. Overall, our model results highlight the crucial role of extreme and compound events in driving severe impacts on pelagic marine ecosystems.

Of particular concern are compound events, which occur when conditions are anomalous for multiple ocean ecosystem drivers (Burger et al., 2022;Gruber et al., 2021;Le Grix et al., 2021, 2022;Zscheischler et al., 2018).The "Blob," a prolonged and extensive marine heatwave that occurred from 2013 to 2015 in the Northeast Pacific, illustrates the potential threat posed by marine compound events on ecosystems.This marine heatwave coincided with anomalously low oxygen, low pH, and large negative anomalies in phytoplankton NPP (Gruber et al., 2021;Le Grix et al., 2021;Mogen et al., 2022;Whitney, 2015;Wyatt et al., 2022), leading to severe impacts on marine life (Cavole et al., 2016), including mortality and reproductive failure of sea birds (Jones et al., 2018;Piatt et al., 2020), mass strandings of sea lions in California and of whales in the western Gulf of Alaska (Cavole et al., 2016), as well as shifts in species distribution toward warm-water species, with repercussions on fisheries (Cavole et al., 2016;Cheung & Frölicher, 2020).Previous research has shown that marine heatwaves and compound extreme events have become more frequent over the past century (Gruber et al., 2021;Oliver et al., 2018) and that this trend is projected to continue as global warming persists (e.g., Burger et al., 2022;Frölicher et al., 2018).If extreme events and compound events regularly induce collapses in animal biomass and community reorganization, the consequences of an increase in their frequency could be catastrophic for ecosystems, fisheries, and human coastal communities.However, the extent to which ocean extreme events and compound events have negative impacts on marine ecosystems remains unclear.
Risks to marine ecosystems arise not only from hazards, such as marine heatwaves or compound events, but also from the interactions between hazards, exposure, and vulnerability (e.g., Bindoff et al., 2019;Magnan et al., 2021).Marine ecosystems may not be exposed nor vulnerable to certain extreme or compound extreme events.For example, Fredston et al. (2023) use a collection of bottom trawl data from Atlantic and Pacific marine ecosystems to show that historical marine heatwaves did not substantially impact the community composition and biomass of these ecosystems.Certain marine species may even benefit from extreme events.Cavole et al. (2016) reported increased recruitment of rockfish in California and northward expansion of tropical and subtropical copepods during the "Blob."These findings highlight the complexity in the relationship between hazards and impacts on marine ecosystems and suggest that compound events with moderate anomalous ocean conditions may also drive or contribute to severe impacts (Zscheischler et al., 2018).To effectively predict and mitigate future impacts on marine ecosystems, a better understanding of the ocean conditions leading to extreme impacts on marine ecosystems is needed.
The limited understanding of the drivers of extreme impacts on marine ecosystems is partly a result of a lack of sufficient observations (Gruber et al., 2021).One approach that circumvents this lack of observations is the use of ensemble climate-impact modeling simulations (Tschumi et al., 2022;van der Wiel et al., 2020).These simulations couple a climate model with an impact model, in our case a global marine fish model.Ensemble climate simulations are produced with the same single climate model under identical external forcing but starting from different initial conditions (Deser et al., 2020;Frölicher et al., 2009).These simulations are then used to force the marine fish model, resulting in a large dataset from which to analyze rare events (e.g., Bevacqua et al., 2023;Cheung et al., 2021;Le Grix et al., 2021;Maher et al., 2021;Poschlod et al., 2020).The ensemble climate-impact modeling simulations are used for two different purposes.The first and most common purpose is forward modeling, which samples oceanic events, such as marine heatwaves, and quantifies their associated impacts on marine ecosystems (Cheung et al., 2021;Cheung & Frölicher, 2020).However, this approach requires prior knowledge of potentially harmful hazards and typically only considers extreme events and ignores other moderate drivers of extreme impacts.Here, we employ a backward approach (e.g., Ben-Ari et al., 2018;van der Wiel et al., 2020;Vogel et al., 2021;Zscheischler, Michalak, et al., 2014;Zscheischler, Reichstein, et al., 2014), which starts by sampling events of extreme impacts and then looks back at the ocean conditions potentially causing these events.This "impact-driven" approach allows for the discovery of unexpected drivers (van der Wiel et al., 2020) of extreme impacts on fish biomass.
The goal of this study is to identify the ocean conditions associated with the most severe impacts on pelagic fish species, especially those associated with extremely low fish biomass.To achieve this, we employ the global marine fish model DBEM, driven by output of a large ensemble simulation of the comprehensive Earth system model GFDL ESM2M.By considering the total biomass of 326 pelagic fish species, we aim to gain a deeper understanding of the drivers of extreme impacts on the entire pelagic community, rather than focusing on positive and negative changes in the biomass of individual species.While this study focuses on a specific climate-fish impact model that is associated with particular assumptions and uncertainties, we aim to obtain generalizable insights that would form the foundation for future studies.

| Ensemble climate-impact modeling
We apply an ensemble climate-impact modeling (van der Wiel et al., 2020;Vogel et al., 2021) approach to identify the environmental drivers that lead to projected low pelagic fish biomass events.The approach consists of three steps: (1) forward modeling; (2) identification of low biomass events; (3) backward assessment of the drivers of low fish biomass events.These three steps are illustrated in Figure 1 and described in detail in the following paragraphs.

| Step: Forward modeling
We use annual mean fish biomass data of 326 pelagic fish species simulated by the species-based Dynamic Bioclimatic Envelope Model (DBEM; Cheung et al., 2008Cheung et al., , 2016) ) (Figure 1).Pelagic fish represent a large proportion of the marine fish biomass relevant to fisheries (Pauly & Palomares, 2020).The DBEM uses a species distribution modeling algorithm (Close et al., 2006) to estimate the initial distribution of exploited species based on their maximum and minimum depth limits, northern and southern latitudinal range limits, habitat type, and known occurrence boundaries (Cheung et al., 2008).These input parameters are mainly provided by two online databases: FishBase (www.fishbase.org) and SeaLifeBase (www.sealifebase.org)and limited to the United Nations Food and Agriculture Organization (FAO) areas.Both the range of FAO areas and the algorithm by Close et al. (2006) constrain pelagic fish distribution in the DBEM, such that pelagic fish species are absent from all regions in white in Figure 2, that is, in the Arctic Ocean and parts of the Southern Ocean.Once the DBEM has determined the species distribution, it then simulates their growth, population dynamics, and net migrations, depending on ocean temperature, oxygen, and advection, as well as a set of species-specific growth parameters.
Movement and dispersal of adults and larvae are modeled through advection-diffusion reaction equations.Environmental preferences are identified for each species by overlaying environmental data from the Earth system model GFDL ESM2M (see below) with maps of the species relative abundance.When environmental conditions deviate from a species environmental preferences, habitat suitability decreases, resulting in a decrease in the species abundance.Biomass is calculated from the population mean body weight and abundance.
The DBEM projects shifts in marine species biomass and distribution under changes in ocean conditions (e.g., Cheung et al., 2013Cheung et al., , 2016) ) including changes in ocean extreme conditions, such as marine heatwaves, that are consistent with alternative species distribution models and empirical evidence, where data exist (e.g., Cheung et al., 2021).The horizontal resolution of the DBEM is 0.5° longitude × 0.5° latitude.The model was spun up for a hundred years using ocean conditions from 1971 to 2000, allowing the population to reach an equilibrium before it was perturbed with environmental changes from 1951 to 2000.
The DBEM was driven by environmental data from a 10-member large ensemble simulation (LES) of an Earth system model, the GFDL ESM2M, covering the time period 1951-2000 (Figure 1).The GFDL ESM2M, developed at NOAA's Geophysical Fluid Dynamics Laboratory (GFDL), is a fully coupled carbon cycle-climate Earth system model (Dunne et al., 2012(Dunne et al., , 2013)), which participated in the Coupled Model Intercomparison Project Phase 5 (CMIP5).It F I G U R E 1 Schematic of the ensemble climate-impact modeling approach used in this study.In a first step ("forward modeling"), changes in fish biomass for 326 pelagic species were simulated over the period 1951-2000 by the Dynamic Bioclimatic Envelope Model (DBEM) forced with environmental data from a 10-member large ensemble simulation of the Earth System Model GFDL ESM2M on a 0.5° longitude × 0.5° latitude grid.In a second step ("Identification of low fish biomass events"), low fish biomass events were identified using the 10th percentile threshold (red areas in the bottom figure).In a third and final step ("backward assessment of the drivers of low fish biomass"), the environmental conditions driving low fish biomass events are assessed.The upper right panel illustrates low fish biomass events sampled at a grid cell in the northern Pacific (at 71° N; 167° W) for one ensemble member, when the pelagic fish biomass is lower than its local 10th percentile computed from 1951 to 2000 over the entire 10-member simulation.
couples an atmospheric circulation model to an oceanic circulation model and includes representations of land, sea ice, and iceberg dynamics, as well as interactive biogeochemistry.The ocean biogeochemical module, TOPAZv2 (Dunne et al., 2013), simulates 30 tracers, including three phytoplankton groups (small and large phytoplankton, diazotrophs) and implicit zooplankton activity.The horizontal resolution of the ocean model MOM4p1 (Griffies, 2012) is nominally 1° with increasing meridional resolution of up to 1/3° toward the equator.The 10-member large ensemble simulation is forced with prescribed historical concentrations of atmospheric CO 2 and non-CO 2 radiative forcing agents over the historical period (Burger et al., 2020(Burger et al., , 2022)).The GFDL ESM2M data were regridded to a horizontal resolution of 0.5° for use in the DBEM.Input GFDL ESM2M data into the DBEM include annual mean horizontal velocities, temperature, dissolved O 2 , and salinity at the surface of the ocean, vertically integrated NPP (sum of small and large phytoplankton and diazotrophs), and sea ice coverage.These DBEM simulations do not take into account changes in acidity (e.g., Tai et al., 2021) or fishing pressure (e.g., Cheung et al., 2018).We force the DBEM with annual mean environmental data for computational efficiency because prolonged anomalies in ocean conditions, like extended high temperatures lasting at least a year, potentially exert the most substantial influence on fisheries (Cheung et al., 2021, Materials and Methods).

| Step: Identification of low fish biomass events
Next, we identify low fish biomass (LFB) events over the period 1951-2000 (Figure 1).We chose that time period as it represents the historical oceanic state and is short enough to not contain too large long-term trends in ocean variables.At each grid cell, we define LFB events as years when the annual mean total biomass of pelagic fish is lower than its local 10th percentile (Figure S1) computed from all ensemble members over 1951-2000.Therefore, these events correspond to 1-in-10-year events.From the 10 realizations of the 50-year period, we thus identify 10 × 5 = 50 LFB events per grid cell.
Defining LFB events with a lower percentile threshold, like the 5th percentile, does not qualitatively change our results.

| Step: Backward assessment of the drivers of low fish biomass events
In a third step, we investigate the drivers of extreme reductions in simulated pelagic fish biomass over the 1951-2000 historical period.
To do this, we perform a backward assessment of the environmental drivers of LFB events (e.g., Ben-Ari et al., 2018;Gagné et al., 2020;van der Wiel et al., 2020;Vogel et al., 2021;Zscheischler, Michalak, et al., 2014;Zscheischler, Reichstein, et al., 2014).Changes in pelagic fish biomass in the DBEM can be driven by changes in depthintegrated net primary production (NPP), surface temperature (T), surface dissolved oxygen levels (O 2 ), surface salinity (S), and sea ice coverage (Ice), that is, the percentage of the grid cell area covered by sea ice.We standardize these ocean variables so as to best compare their contributions to driving LFB events, by removing the mean and dividing by the standard deviation computed over 1951-2000.
To identify drivers of LFB events, we employ a LASSO (least absolute shrinkage and selection operator) logistic regression (Tibshirani, 1996;Vogel et al., 2021).This statistical model allows for classifying years into LFB events and non-LFB events depending on ocean conditions.We consider ocean conditions up to 2 years before the event, to account for possible lagged effects.The probability of an LFB event to occur is given as: with where X i [NPP, T, O 2 , S, and Ice] at the year of the event, 1 year prior to the event (indicated with −1y), and 2 years prior to the event (indicated with −2y).An LFB event is predicted when P(LFB) > 0.5.
The regression coefficient i accounts for the link between an ocean variable X i and the probability of LFB.A positive i signifies that an increase in X i raises the probability of LFB, and vice versa.However, high correlation between NPP, T, O 2 , S, and Ice implies a high variability of the coefficients i (Vogel et al., 2021).
For example, T and O 2 are often negatively correlated in the surface ocean, where the main driver of oxygen changes is oxygen solubility, which decreases with increasing temperature (Garcia & (1) Area under the precisionrecall curve (AUC PR) of the LASSO logistic regression.Model performance is satisfactory over areas where AUC PR > 0.1.White areas correspond to regions without fish biomass data.Gordon, 1992).The information brought by a high absolute value of T and a low absolute value of O 2 could alternatively be conveyed by a low absolute value of T and a high absolute value of O 2 .To address the high correlation between variables, we specifically employ a LASSO logistic regression (Tibshirani, 1996), which prevents high variability in the coefficients by applying a penalty term on the norm of the coefficients ‖ ‖.The regression coefficients are determined by minimizing a cost function Cost ( ) with regularization on ‖ ‖: The parameter controls the strength of the regularization.
Through fivefold cross-validation, we obtained optimal , the value of associated with the highest mean cross-validated performance of the model.We then selected = SE , the value of that gives the most regularized model (i.e., the highest λ) such that its crossvalidated performance is within one standard deviation of that of optimal (Friedman et al., 2010;Krstajic et al., 2014).The penalty term ‖ ‖ tends to produce some regression coefficients i that are exactly 0. The LASSO logistic regression therefore performs an automatic selection of the variables that are statistically linked to LFB.
As a result, certain grid cells might have just one predictor, like NPP, while others might have up to 15 predictors.
Model performance is assessed by the area under the precisionrecall curve (AUC PR), a metric commonly used to summarize the performance of a binary classification model when the sets are unbalanced (here, 10% of years are LFB events, 90% are not; Cook & Ramadas, 2020;Saito & Rehmsmeier, 2015).More than 99.99% of the ocean surface area considered by the model, excluding regions without modeled fish biomass (white areas in Figure 2) has an AUC PR value greater than 0.1 (Figure 2), the rate of LFB events in the time series (Saito & Rehmsmeier, 2015), which indicates that the logistic regression performs well at predicting LFB events.Note that we allowed for the inclusion of lagged effects.Namely, variables Allowing for lagged effects increases model performance by 35% compared with a model based on concurrent predictors only.

| Categorizing the drivers of low fish biomass events
We use the coefficients from the LASSO logistic regression to assess whether extreme, compound, or compound and extreme drivers are necessary conditions for LFB events.An LFB is predicted when P(LFB) > 0.5, which is equivalent to y > 0 in Equation ( 2).For a given location and given the regression coefficients i , we test whether y can be positive with all predictors between their 10th and 90th percentile.If not, moderate events cannot drive LFB events, and extreme events are necessary to drive LFB events.Similarly, we test whether drivers can be univariate by testing whether y can be positive with only one predictor being nonzero (and the rest are zero).If not, compounding drivers are a necessary condition for LFB events.

| RE SULTS
Low fish biomass events are associated with shifts in ocean conditions, as demonstrated by the red and blue distributions presented in Figure 3. On a global scale and on average, NPP is decreased by 0.8 standard deviation (SD) during LFB events compared with normal conditions.Temperature is increased by 0.5 SD, dissolved oxygen decreased by 0.5 SD, and salinity is 0.1 SD lower than usual.Furthermore, sea ice coverage variability is also enhanced during LFB events.In the following subsections, we further discuss these changes in ocean conditions and show how they can be used to predict LFB events.

| Drivers of low fish biomass events
Using LASSO logistic regression, the probability of an LFB event to occur is modeled as a function of NPP, T, O 2 , S, and Ice conditions up to 2 years prior to the event.Figure 4a presents the global mean regression coefficients associated with each predictor.A positive coefficient indicates that any increase in the associated predictor increases the probability of an LFB event, while a negative coefficient indicates that any decrease in the associated predictor increases the probability.We first describe the regression coefficients associated with NPP, T, O 2 , S, and Ice conditions in the year of the event.
Globally, the LASSO logistic regression coefficients indicate that LFB events are more likely to occur when NPP, O 2 , and S are anomalously low, and when T and Ice are anomalously high (Figure 4a).This is consistent with the shift in these ocean variables' distribution during LFB events (Figure 3).NPP is selected over 88% of the global ocean (Figure 4b), due to its negative impacts on the biomass of fish species in the DBEM (Cheung et al., 2011).On the contrary, T, O 2 , and S are selected over 47%, 42%, and 29% of the global ocean, respectively (Figure 4b).The effect of changes in T on fish biomass varies depending on the species' temperature preference, which may result in poor prediction of LFB events when considering the total biomass of all pelagic fish species.Changes in salinity are usually too subtle in the open ocean to directly affect osmotic processes in fish and ultimately limit fish biomass.However, salinity changes may be correlated with changes in other ocean conditions that directly influence fish biomass and thus serve as a proxy for LFB events in certain regions.Ice is only selected as a predictor in and around sea-ice-covered regions (5% of the global ocean, Figure 4b).
There is a substantial degree of spatial variability, as indicated by the standard deviation of each regression coefficient over the area over which it is selected (Figure 4a).The regression coefficient for NPP (Figure 5a) is negative over most of the global ocean, which means that low NPP is related to LFB.The more negative the regression coefficient, the more a reduction in NPP increases the likelihood of LFB.For example, in the northeastern Pacific, where the regression coefficient for NPP is particularly negative, even a relatively weak reduction in NPP might be sufficient to drive LFB events.In contrast, in the center of the northern Pacific, a stronger reduction in NPP might be necessary, as indicated by the less negative regression coefficient.
LFB events are associated with high temperatures over most of the area over which T is a predictor (positive coefficient in Figure 5b), yet they are also associated with low temperatures over a few regions in the northern high latitudes and in the northern part of the Southern Ocean.The regression coefficient for O 2 varies in sign over the global ocean (Figure 5c).Over most of the low-to mid-latitudes, LFB events are associated with low O 2 .Low oxygen levels limit metabolism, growth performance and body size, and therefore limit total population biomass (Cheung et al., 2011(Cheung et al., , 2013;;Clarke et al., 2021;Pauly & Cheung, 2017), potentially driving LFB events.In isolated grid cells in the high northern latitudes and in the northern part of the Southern Ocean, LFB events are, however, associated with high O 2 .There, high O 2 is usually driven by increased solubility associated with lower T (Figure S3), a potential driver of LFB for which high O 2 would simply be a proxy.Low salinity is an indicator of LFB especially in the equatorial Pacific, whereas high salinity is an indicator of LFB in the western South Indian Ocean and the North Atlantic (Figure 5d).LFB events are also favored by high Ice in the northern high latitudes and low Ice in the southern high latitudes (Figure 5e).This hemispheric divergence is explained by the species composition of the pelagic ecosystem in each region.In the Arctic, although larger sea ice coverage provides more favorable conditions for polar species, it reduces the suitable habitat for nonpolar, nonendemic species (Cheung et al., 2008).Reduced biomass of nonpolar species during high Ice events reduces the overall pelagic biomass and contributes to driving LFB events in the Arctic.In contrast, nonpolar species are less connected with the Antarctic, where marine species are highly endemic and therefore well adapted to sea ice (Eastman, 2005).Reduced biomass of polar species during low Ice events contributes to driving LFB events around Antarctica in the DBEM.

| The role of lagged effects in driving low fish biomass events
Anomalous ocean conditions may have lagged effects on fish biomass, which potentially drive LFB events a few years later.To account for these lagged effects, ocean conditions 1-2 years prior to an LFB event were also selected as predictors in our analysis.
We found that NPP from 1 and 2 years prior to an LFB event was selected as a predictor over about 70% of the ocean area considered by the model (Figure 4b).The associated regression coefficients are still relatively high (on the order of −0.3; Figure 4a).This lagged effect might be explained by NPP having large negative impacts on fish population biomass (Chassot et al., 2010), leading to a reduction in reproductive capacity that would take multiple years to recover particularly for longer-lived, later-matured species (i.e., with lower intrinsic population growth rate).Negative impacts from low NPP propagate over time through the population, potentially driving a decline in overall fish biomass.The lagged effect is most pronounced in the low latitudes (Figure S4), where time variability in NPP is also especially low in the ESM2M (Le Grix et al., 2022).Low NPP may therefore be associated with low NPP in the following years and thus indirectly with LFB in the following years.A time lag of 1 year or 2 years for T and O 2 is rarely selected as a predictor (over less than 21% of the ocean area, Figure 4b).Moreover, the regression coefficients associated with lagged T and O 2 are of much lower absolute value compared with the regression coefficients associated with concurrent T and O 2 (Figure 4a).T and O 2 appear better suited F I G U R E 5 Regression coefficient associated with net primary production (a), surface temperature (b), dissolved oxygen (c), salinity (d), and sea ice coverage (e) conditions in the year of the event in the logistic regression.Black areas correspond to regions without fish biomass data, in contrast to white areas, which correspond to regions where a given ocean variable was not selected as a predictor in the logistic regression.
at predicting concurrent LFB events than future LFB events.Compared with concurrent salinity, one-year and two-year lagged S was selected as a predictor over a similar area (30% of the ocean area considered by the model).This suggests persistent low salinity conditions, which indicate negative impacts on fish biomass over time.
In contrast, one-year and two-year lagged Ice is never selected as a predictor, potentially reflecting high interannual variability in sea ice coverage.Maps of the regression coefficient associated with oneyear and two-year lagged predictors are shown in Figure S4.
Our results highlight the important role of lagged effects in driving LFB events.They might include "direct" lagged effects, by which anomalous ocean conditions drive LFB on their own after 1 or 2 years, as well as "indirect" lagged effects.Anomalous ocean conditions may indirectly drive LFB by playing the role of a preconditioning event or by temporally compounding with other events (Zscheischler et al., 2020).Preconditioning events correspond to anomalous ocean conditions rendering an ecosystem more vulnerable to subsequent ocean events, which might end up driving LFB.
Temporally compounding events correspond to a repetition of the same kind of ocean event, like a marine heatwave, whose impacts accumulate over time and end up driving LFB.

| Most influential predictor of low fish biomass events
We identified at each grid cell the most influential predictor of LFB events, that is, the predictor associated with the highest absolute regression coefficient (Figure 6).NPP is the main predictor over most (52%) of the ocean area considered by the model.−1-year and −2-year lagged NPP are the main predictors over 16% of the ocean, mostly located in the subtropical gyres, where NPP variability is low.T is the main predictor in the eastern equatorial Pacific, in the equatorial Atlantic and locally over the Gulf Stream region and in the high latitudes.O 2 is the main predictor over a few regions, for example, locally in the eastern equatorial Pacific and Atlantic.S is the main predictor in the western equatorial Pacific.Lastly, Ice is the main predictor over certain regions in the high latitudes, such as in the Arctic Ocean.
Overall, our results highlight the dominance of NPP, as a useful indicator of potential extreme declines in pelagic fish biomass over most of the ocean area (Figure 6).This is consistent with previous studies on the link between long-term changes in NPP and fish production (e.g., Blanchard et al., 2012;Chassot et al., 2010;Lotze et al., 2019;Stock et al., 2017), although we look here at annual changes.

| Categorizing the drivers of low fish biomass events into moderate, extreme, univariate, and compound drivers
We use LASSO logistic regression with lagged effects to classify the drivers of LFB events into four categories: univariate, compound, moderate, and extreme.Our methodology utilizes the regression coefficients from the previously built logistic regression to determine whether an LFB event can be driven by univariate or moderate anomalies in ocean conditions or if compound or extreme anomalies are necessary (see Section 2).
Our analysis shows that moderate oceanic events can drive LFB events over only 30% of the ocean area considered by the model (in gray in Figure 7a), primarily in the low latitudes, particularly in F I G U R E 6 Most influential predictor of LFB events, defined as the predictor with highest absolute regression coefficient in the logistic regression.These regions have a high number of predictors in the logistic regression, and moderate anomalies in multiple variables may compound to drive LFB events.Extreme oceanic events are necessary to experience LFB events over 70% of the ocean area, primarily in the midto-high latitudes (yellow area in Figure 7a).This is further supported by the co-occurrence of LFB events with extreme conditions, such as extremely high SST, low O 2 , and low NPP in the mid-latitudes, as well as with extremely high Ice and low Ice in the northern and southern high latitudes, respectively (Figure S5).
Next, we analyze whether univariate events alone can drive LFB events or if compound events are necessary.Our results show that univariate events can drive LFB events over 22% of the ocean area, primarily in the eastern equatorial Pacific (in gray in Figure 7b), where high absolute regression coefficient for O 2 and T suggests that even univariate O 2 and T events may be sufficient to drive LFB events.
Over the remaining 78% of the ocean area (in blue in Figure 7b), compound events are necessary to experience LFB events.In the Appendix, we further differentiate between temporally compounding events, where one single ocean variable is anomalous over multiple years, and multivariate compound events, where multiple ocean variables are anomalous (Figure S6).Temporally compounding events are necessary for the occurrence of LFB events over 36% of the ocean area (in light blue in Figure S6), mostly located in the subtropics, where NPP, −1 year NPP, and −2 year NPP are LFB predictors.There, low NPP could persist over multiple years and drive LFB.In contrast, multivariate compound events are necessary for LFB events over the remaining 42% of the ocean area (in dark blue in Figure S6), where multiple events such as high T, low O 2 , and low NPP events may compound to drive LFB.
Lastly, we overlapped Figure 7a,b to identify regions where the LFB drivers must be both an extreme and a compound event (Figure 7c).The fraction of ocean area where extreme and compound events are necessary to experience LFB events overlap in the mid-to-high latitudes (61% of the ocean, in green in Figure 7c).
There, a compound extreme event with an extreme in at least one variable is required to drive extreme impacts on pelagic fish.
Note that although compound extreme events are necessary to experience LFB events over 61% of the ocean area (in green in Figure 7c), this area comprises only 37% of the total pelagic fish biomass (also in green in Figure S7c), as fish species are heterogeneously distributed over the global ocean (Figure S1a).Over the remaining 13% of the ocean area (in gray in Figure 7c), LFB events can be driven by a univariate event as long as it is extreme or, interchangeably, by a moderate event as long as it is a compound moderate event (i.e., a combination of moderate anomalies in multiple ocean variables).A univariate and moderate anomaly in one ocean variable is not sufficient to drive LFB anywhere (in black on the legend of Figure 7c).In the appendix, we converted the fraction of ocean area occupied by each driver category (Figure 7) into the corresponding portion of the global mean fish biomass within each category (Figure S7), addressing heterogeneous fish biomass over the global ocean.

| D ISCUSS I ON AND CON CLUS I ON
We investigated the drivers of extremely low pelagic fish biomass (LFB) events by applying a LASSO logistic regression to output of a global marine fish model.We found that low net primary productivity is the most influential predictor of LFB events over the majority (68%) of the surface ocean considered by the model.The prediction of such LFB events is substantially improved by considering net primary production 1-2 years before the event.We also found that a moderate and univariate anomaly in one ocean variable is not sufficient to drive LFB events and that over 61% of the surface ocean, anomalies in multiple ocean variables-which must be extreme in at least one variable-are required to cause LFB events.Our findings highlight the key role of extreme and compound events for severe impacts on pelagic fish biomass.
This study takes the original approach of investigating the drivers of extreme impacts on pelagic fish biomass using a backward assessment method.Contrary to previous studies that apply a forward approach to investigate the impacts of extreme hazards, such as marine heatwaves, on ecosystems (e.g., Cheung et al., 2021;Cheung & Frölicher, 2020;Smale et al., 2019;Smith et al., 2023), the focus here is rather on extreme impacts on fish and what drives such impacts.
Such backward approaches have become more common in recent years, for instance to identify drivers of extremely low vegetation activity and carbon uptake (Zscheischler et al., 2013;Zscheischler, Mahecha, et al., 2014), floods (Jiang et al., 2022), and crop failure events (Vogel et al., 2021).The backward approach allows for the identification of potentially unexpected drivers (van der Wiel et al., 2020), such as low NPP, and for the distinction between univariate, compound, moderate, and extreme drivers.The dominant role of low NPP, and not of high temperature (e.g., Cheung & Frölicher, 2020), in driving extremely low fish biomass, was particularly unexpected.
Our findings may explain why some recent studies did not find substantial impacts of marine heatwaves on ecosystems, for example, on the biomass and community composition of coastal Atlantic and Pacific demersal fish communities over the last three decades (Fredston et al., 2023).In contrast, our results are supported by previous studies arguing that long-term changes in fish production strongly mirror long-term changes in NPP (e.g., Blanchard et al., 2012;Chassot et al., 2010;Lotze et al., 2019;Stock et al., 2017).Marine heatwaves possibly lead to negative or more negative impacts when compounded with other stressors, such as low NPP.A multivariate approach is recommended to best apprehend changes in fish biomass.Additionally, we found that combinations of moderate ocean anomalies might be sufficient to drive extreme impacts on pelagic fish biomass over 22% of the ocean area considered by the model.Therefore, monitoring multiple ocean ecosystem drivers in addition to marine heatwaves (e.g., Jacox et al., 2022) may help to improve predictions of high impacts on the pelagic fish ecosystem.
An important aspect of this study is that it considers the perspective of a very large assemblage of pelagic fish species, as opposed to focusing on individual species.Note that our analysis excludes demersal fishes and benthic invertebrates.Each pelagic species has its own unique response to anomalous ocean conditions, based on its natural habitat and range of tolerance (Pörtner, 2002;Pörtner & Peck, 2010).For instance, some species that are already at the upper limit of their thermal tolerance might be negatively affected by elevated temperatures, while others may benefit from the same warming event.However, the overall response of the pelagic fish ecosystem to changes in ocean conditions is not yet well understood.This study offers an initial evaluation of the drivers of extreme declines in the biomass of pelagic fishes.It forms the foundation for future studies, which could then focus on specific oceanic events, known to impact pelagic ecosystems.Such studies could be conducted at a finer temporal scale, so as to also address the role of event intensity and duration in driving impacts.
A compound event combines individual ocean events, whose effects may act synergistically to drive LFB events (Boyd & Brown, 2015;Gruber et al., 2021).We employed a logistic regression to analyze these effects, assuming an additive relationship between them.However, it is important to note that these effects may not be linear, as demonstrated by previous research (Zscheischler & Seneviratne, 2017).For example, the combined effect of increased temperature and reduced net primary production on fish may be different from the sum of their isolated effects.To further investigate this, we conducted a sensitivity analysis (not shown) that considered mixed effects in the logistic regression.Specifically, we allowed the product of two ocean variables to be selected as a predictor of LFB events under the condition that it improved the prediction.However, in our study, accounting for mixed effects did not improve the prediction of LFB events.In contrast, it reduced the model performance by 14%, possibly due to difficulties in selecting the best predictors among many variables.
The robustness of the results presented here depends on the fidelity of the DBEM in simulating the ocean ecosystem drivers of changes in pelagic fish biomass.The DBEM has been shown to reproduce observed shifts in the biomass and distribution of hundreds of fish species under moderate (e.g., Cheung et al., 2013Cheung et al., , 2016) ) and extreme (Cheung et al., 2021) changes in ocean conditions.However, as with any model, the DBEM is an incomplete representation of the pelagic ecosystems (Cheung et al., 2011) and certain limitations require further discussion.
One limitation of the DBEM is that it does not account for vertical heterogeneity in ocean ecosystem drivers and assumes surface ocean conditions to be the primary drivers of pelagic fish biomass.In reality, pelagic fishes are also distributed at subsurface, where changes in ocean conditions do not necessarily mirror changes in surface ocean conditions.This is particularly true in not well-ventilated regions, such as oxygen minimum zones, where oxygen changes are decoupled from changes in surface oxygen.Further studies are needed to better assess the role of subsurface oxygen anomalies in driving extreme declines in pelagic fish biomass (e.g., Morée et al., 2023).Another limitation of the DBEM is that it integrates the effect of recruitment and mortality through an intrinsic population growth model; therefore, we cannot disentangle the respective contributions of a change in recruitment or mortality to reducing fish biomass during LFB events.Additionally, the DBEM simulations that we used in this study do not consider the effects of ocean acidification or fishing pressure on fish biomass.Ocean acidification can cause physiological stress and perturb fish olfactory ability to detect suitable habitat, potentially leading to population declines (Branch et al., 2013;Cheung et al., 2011;Melzner et al., 2009;Munday et al., 2009;Tai et al., 2021).These effects may be particularly pronounced when combined with other ocean stressors such as low oxygen (Gobler & Baumann, 2016) or high temperature (Burger et al., 2022;Cornwall et al., 2021).Fishing pressure also impacts fish biomass (e.g., Watson et al., 2013) and may be a primary reason for the recent marine biodiversity decline (Jaureguiberry et al., 2022).
However, fisheries target specific species, and each species reacts differently not only to the fishing pressure but also to environmental stressors.Therefore, one cannot disentangle the respective effects of fishing and the environment on the whole pelagic ecosystem.The DBEM offers the possibility to remove the fishing pressure and concentrate on the environmental drivers of LFB events.Once the environmental drivers are understood, further studies could address the additional role of fisheries in driving extreme changes in fish biomass (Cheung et al., 2011(Cheung et al., , 2022;;Tai et al., 2018).
The conclusions of this study are not only dependent on the global fish model used but also on the ocean ecosystem drivers as simulated by the GFDL ESM2M.For example, the dependencies between ocean variables, such as the correlation between net primary production and temperature, may differ between the GFDL ESM2M and other Earth system models.Dependencies between ocean variables are reflected in the logistic regression coefficients and may impact our results.In particular, the negative correlation between net primary production and temperature in GFDL ESM2M is overestimated in the tropics and strongly underestimated around Antarctica when compared to observation-based data (Le Grix et al., 2022).
Therefore, the role of compound high temperature and low net primary productivity events in driving low fish biomass events may be overestimated in the tropics and underestimated around Antarctica.
In addition, the GFDL ESM2M has been shown to underestimate NPP in sea-ice-covered regions (Dunne et al., 2013).Sea ice limits incoming solar radiation and therefore phytoplankton growth in the model, although previous studies suggest sustained NPP within the sea ice cover by sea-ice-bound algae, which support a thriving polar ecosystem (Eicken, 1992;Koch et al., 2023).To further validate the robustness of our results, different Earth system models with potentially divergent dependencies between ocean variables should be used, following similar sensitivity experiments as over land (Tschumi et al., 2023).New simulations from global fish models forced by different Earth system models, that have now become available under the Inter-Sectoral Impact Model Intercomparion Project (ISIMIP) framework (Tittensor et al., 2021), might be used for an intermodel comparison study in the future.In this study, the selection of the GFDL ESM2M was motivated by the availability of a large ensemble simulation (LES), which provides the necessary large dataset from which to study rare extreme events.At present, there are only a few LES available with fully coupled Earth system models that simulate ocean ecosystem drivers (e.g., Deser et al., 2020;Maher et al., 2019;Rodgers et al., 2015Rodgers et al., , 2021)).Our analysis could, for example, be replicated using the 100-member CESM2-LES (Rodgers et al., 2021).
The drivers of LFB events identified in this study are affected by climate change, as the ocean is warming and losing oxygen, sea ice is retreating and primary production is responding to changes in temperature, nutrients, light levels, and grazing pressure (Kwiatkowski et al., 2020).Changes in these environmental stressors are bound to affect marine ecosystems (e.g., Cheung et al., 2013;Doney et al., 2012;Lotze et al., 2019), and therefore likely also the occurrence of LFB events.For example, over regions where high temperatures drive LFB, for example, in the eastern equatorial Pacific and equatorial Atlantic, one would expect an increase in LFB events under global warming.
Similarly, a change in NPP's mean state (Steinacher et al., 2010;Stock et al., 2017) and variability would alter the likelihood of LFB events given the dominant role of low NPP in driving LFB over most of the global ocean.In conjunction with short-term forecasts and longerterm projections of ocean conditions, our findings help anticipate potential declines in pelagic fish biomass.Early prediction allows mitigating measures to be quickly implemented to moderate impacts on ecosystems and the societies relying upon them.
This study provides insights into potential environmental drivers of high impacts on pelagic fish biomass.We highlight the key role played by univariate extreme events, compound moderate events, and compound extreme events in driving severe impacts on fish.Our results motivate further work on ocean extreme and compound events, as well as the monitoring of multiple ocean ecosystem drivers including net primary productivity as a means to predict impacts on pelagic fish.
1 year (−1 year NPP, −1 year T, −1 year O 2 , −1 year S, −1 year Ice) and 2 years (−2 years NPP, −2 years T, −2 years O 2 , −2 years S, −2 years Ice) before the event can also be selected as predictors.The LASSO logistic regression only selects as predictors the variables that are statistically linked to LFB.It may therefore select, for example, −1 year NPP and/or −2 years NPP but exclude −1 year O 2 and/or −2 years O 2 .
Global distribution of the annual mean standardized (a) net primary production (NPP), surface (b) temperature (T), (c) dissolved oxygen (O 2 ), (d) salinity (S), and (e) sea ice coverage (Ice) over all years (in blue) and low fish biomass years (in red) over 1951-2000.On the yaxis, density accounts for the probability density function, whose area under the curve is equal to 1. On the x-axis, a shift of the distribution toward −1 corresponds to a reduction in annual mean values of one standard deviation.F I G U R E 4 (a) Global mean logistic regression coefficient associated with the net primary production (NPP), surface temperature (T), dissolved oxygen (O 2 ), salinity (S), sea ice coverage (Ice), −1 year NPP (1 year prior to the event), −1 year T, −1 year O 2 , −1 year S, −1 year Ice, −2 years NPP, −2 years T, −2 years O 2 , −2 years S, and −2 years Ice.Error bars indicate the standard deviation of each regression coefficient across space.An error bar crossing the x-axis signifies that low fish biomass (LFB) events are associated with a positive or negative anomaly of a given predictor, depending on the region.(b) Fraction of the ocean area considered by the model (%) over which each variable is selected as a predictor of LFB events in the regression.

F
I G U R E 7 Spatial distribution of each category of LFB drivers.The legend indicates the fraction of ocean area considered by the model occupied by each category.(a) Extreme events are necessary to experience LFB events.(b) Compound events are necessary to experience LFB events.Both upper panels are superimposed to create the bottom panel (c), which shows where compound extreme events are necessary to experience LFB events.the central equatorial Pacific, and in the California Current System.