Including ﬁlter-feeding gelatinous macrozooplankton in a global marine biogeochemical model: model-data comparison and impact on the ocean carbon cycle

Filter-feeding gelatinous macrozooplankton (FFGM), namely salps, pyrosomes and doliolids are increasingly recognized as an essential component of the marine ecosystem. Unlike crustacean zooplankton (eg., copepods) which feed on prey that is an order of magnitude smaller, ﬁlter-feeding allows FFGM access to a wider range of organisms, with predator over prey ratios as high as 100 000:1. In addition, most FFGM produce carcasses and/or fecal pellets that sink 10 times faster than those of copepods. This implies a rapid and eﬃcient export of organic matter to depth. Even if these organisms represent < 5% of the overall planktonic biomass, the induced organic matter ﬂux could be substantial. Here we present a ﬁrst estimate of the inﬂuence of FFGM organisms on the export of particulate organic matter to the deep ocean based on a marine biogeochemical earth system model: NEMO-PISCES. In this new version of PISCES, two processes characterize FFGM: the preference for small organisms due to ﬁlter feeding, and the rapid sinking of carcasses and fecal pellets. To evaluate our modeled FFGM distribution, we compiled FFGM abundances observations into a monthly biomass climatology using a taxon-speciﬁc conversion. FFGM contribute strongly to carbon export at depth (0.4 Pg C / yr at 1000m), particularly in low-productivity region (up to 40% of POC export at 1000m) where they dominate macrozooplankton by a factor of 2. This export increases in importance with depth, with a simulated transfer eﬃciency close to one.

Overall, most studies to date have focused on regional scales. Recently Luo et al. (2020) have estimated the contribution to 60 the global carbon cycle of three categories of gelatinous zooplankton: ctenophores, cnidarians and pelagic tunicates. Using a data-driven carbon cycle model, they found that pelagic tunicates contribute three quarters of the particulate organic carbon (POC) flux induced by gelatinous zooplankton or one quarter of the total POC exported at 100 m. A more recent study by the same team (Luo et al., 2022) revised this estimate to 0.57 Pg C yr −1 , representing 9% of total export at 100 m, by explicitly representing FFGM in the COBALT-v2 biogeochemical model (FFGM refer to large pelagic tunicates in their study). 65 Marine biogeochemical models have repeatedly shown their usefulness in understanding marine processes on a global scale: in particular for the role of plankton in ecosystem processes (e.g., Sailley et al., 2013;Le Quéré et al., 2016;Kearney et al., 2021) and biogeochemical fluxes (e.g., Buitenhuis et al., 2006;Kwiatkowski et al., 2018;Aumont et al., 2018). Their complexity has been increased by the addition of multiple limiting nutrients and multiple functional groups or size classes of phytoplankton and zooplankton (e.g., Le Quéré et al., 2005;Follows et al., 2007;Ward et al., 2012;Aumont et al., 2015). 70 Notably, Plankton Functional Type (PFT) models have been introduced as a way of grouping organisms that keeps the overall biological complexity at a manageable level (Moore et al., 2001;Gregg et al., 2003;Le Quéré et al., 2005). Wright et al. (2021) showed that the introduction of a jellyfish PFT (cnidarians only) into the PLANKTOM model has a large direct influence on the biomass distribution of the crustacean macrozooplankton PFT and indirectly influences the biomass distributions of protozooplankton and mesozooplankton through a trophic cascade. This influence could be explained by the specific diet of 75 jellyfish that differs from other zooplankton PFTs. Similarly, the inclusion of FFGM as a new PFT in a PFT-based model has been recently achieved by Luo et al. (2022). In their study, they introduced two tunicates groups into the COBALT-v2 model: a large salps/doliolids (FFGM) and a small appendicularian, and they estimated their impact on the surface carbon cycle, but they did not consider the impacts on the deeper carbon cycle.
Here, we use the PISCES-v2 model (Aumont et al., 2015) which is the standard marine biogeochemistry component of 80 NEMO (Nucleus for European Modelling of the Ocean). In this study, a new version of PISCES was developed (PISCES-FFGM) in which two new PFTs were added: a generic macrozooplankton (GM) based on an allometric scaling of the existing mesozooplankton and a filter-feeding gelatinous macrozooplankton (FFGM). Two processes characterize FFGM in this version of the model: access to a wide range of preys through filter feeding and rapid sinking of carcasses and fecal pellets. We first examine how the model succeeds in reproducing the surface distribution of FFGM by providing a new compilation of 85 abundance observations converted to carbon biomass via taxon-specific conversion functions to make this assessment. Second, because the modeling study by Luo et al. (2022) focused on the impact of FFGM on surface processes, we investigated whether our modeling framework and formulations produce results consistent with theirs. Our study provides also some new insights: 1) we explore the FFGM-specific spatial patterns of organic matter production, export and particles composition in the top 100 m; 2) we investigate the impacts of FFGM on the export of particulate organic carbon to the deep ocean via an explicit 90 representation of fast-sinking fecal pellets and carcasses. , two phytoplankton groups (Diatoms and Nanophytoplankton, denoted D and N ), two zooplankton size classes (Micro-and Mesozooplankton, denoted Z and M ) and an explicit representation of particulate and dissolved organic matter, reaching a total of 24 prognostic variables (tracers). A full description of the model is provided in Aumont et al. (2015). organisms (hereafter referred to as GM, see fig. 1) and the other to salp-like filter-feeding gelatinous macrozooplankton organisms (hereafter referred to as FFGM, see fig. 1). As with micro-and mesozooplankton in the standard version of PISCES, the C:N:P stoichiometric composition of the two macrozooplankton groups is assumed to be constant. In addition to their carbon biomass, two additional tracers were introduced into the model for each macrozooplankton group corresponding to fecal pellets and carcasses in carbon units, respectively (GM Carcasses,GM Fecal Pellets,FFGM Carcasses and FFGM Fecal Pellets,105 see fig. 1). Because both macrozooplankton groups have a constant Fe:C stoichiometry and feed on phytoplankton that have a flexible Fe:C stoichiometry (Eq. 16 to 20 in Aumont et al. (2015)), two compartments representing the iron content of the fecal pellets of the two macrozooplankton groups were added. Figure 1 summarizes the tracers and interactions newly introduced into PISCES for this study (referred to as PISCES-FFGM hereafter).
In total, the tracers considered for particulate and dissolved organic matter are (organic particles in fig. 1): sP OC which 110 refers to the carbon content of small organic particles, bP OC which refers to the carbon content of large organic particles, DOC which refers to dissolved organic carbon, Ca F F GM which refers to the carbon content of FFGM carcasses, F p F F GM which refers to the carbon content of FFGM fecal pellets, Ca GM which refers to the carbon content of GM carcasses and F p GM which refers to the carbon content of GM fecal pellets.

Macrozooplankton (FFGM and GM) dynamics
We first present the generic equation describing the dynamics of the two groups of macrozooplankton, and then focus on the modeling choices differentiating these two groups. All symbols and definitions are summarized in Table 1.
The temporal evolution of the two compartments of macrozooplankton is governed by the following equation: This equation is similar to the one used for micro-and mesozooplankton in PISCES-v2 (Aumont et al., 2015). In this equation, X is the considered macrozooplankton biomass (GM or F F GM ), and the three terms on the right-hand side represent growth, quadratic and linear mortalities. e X is the growth efficiency. It includes a dependence on food quality as presented in PISCES-v2 (Eq. 27a and 27b in Aumont et al. (2015)). Quadratic mortality is divided into mortality due to predation by 125 unresolved higher trophic levels (with a rate m X ) and mortality due to disease (with a rate m X c ). All terms in this equation were given the same temperature sensitivity f X (T ) using a Q10 of 2.14 (Eq. 25a and 25b in Aumont et al. (2015)), as for mesozooplankton in PISCES-v2 and according to Buitenhuis et al. (2006). Growth rate and quadratic mortality are reduced and linear mortality is enhanced at very low oxygen levels, as we assume that macrozooplankton are not able to cope with anoxic waters (∆(O 2 ) varies between 0 in fully oxic conditions and 1 in fully anoxic conditions, see Eq. 57 in Aumont et al. The difference between the two macrozooplankton groups lies in the description of the term G X , i.e. the ingested matter. A full description of the equations describing G X is provided in the appendix Text A3 (Eq. A1 to Eq. A12). Below, we present the two different choices of feeding representation that differentiate the dynamics of the two macrozooplankton groups, GM and FFGM.

135
GM, namely generic macrozooplankton, is intended to represent crustacean macrozooplankton, such as euphausids or large copepods. The parameterization is similar to that of mesozooplankton (Eq. 28 to 31 in Aumont et al. (2015)). Therefore, in addition to conventional suspension feeding based on a Michaelis-Menten parameterization with no switching and a threshold (Eq. A1, A2 and A3), flux-feeding is also represented (Eq. A4) as has been frequently observed for both meso-and macrozooplankton (Jackson, 1993;Stukel et al., 2019). GM can flux-feed on small and large particles as well as on carcasses and fecal 140 pellets produced by both GM and FFGM (Eq. A6). We assume that the proportion of flux-feeders is proportional to the ratio of potential food available for flux feeding to total available potential food (Eq. A7 and A8). Suspension feeding is supposed to be controlled solely by prey size, which is assumed to be about 1 to 2 orders of magnitude smaller than that of their predators (Fenchel, 1988;Hansen et al., 1994). Thus, GM preferentially feed on mesozooplankton, but also, to a lesser extent on microzooplankton, large phytoplankton and small particles (Eq. A5 and A10, Fig. 1).

145
FFGM represent the large pelagic tunicates (i.e. salps, pyrosomes and doliolids but not appendicularians). Pelagic tunicates are all highly efficient filter feeders and thus have access to a wide range of prey sizes, from bacteria to mesozooplankton (Acuña, 2001;Sutherland et al., 2010;Bernard et al., 2012;Ambler et al., 2013). There is no strong evidence that FFGM feed on mesozooplankton in the literature. Therefore, we assume in our model that FFGM are solely suspension feeders (i.e. with concentration-dependent grazing based on a Michaelis-Menten parameterization with no switching and a threshold, see Eq.

Carcasses and fecal pellet dynamics:
Carcasses Ca F F GM and Ca GM are produced as a result of non predatory quadratic and linear mortalities of GM and FFGM, respectively. The F p F F GM and F p GM are produced as a fixed fraction of the total food ingested by the two macrozooplankton 155 groups. Remineralization of fecal pellets and carcasses by bacteria is modeled using the same temperature-dependent specific degradation rate with a Q 10 of 1.9, identical to that used for small and large particles. In addition to remineralization, carcasses and fecal pellets undergo flux feeding by GM as explained in the previous subsection. Note that parasitism is not considered in this study because it is too poorly documented, but that it could represent an important source of carcasses (Lavaniegos and Ohman, 1998;Phleger et al., 2000;Hereu et al., 2020). The sinking speeds of these particle pools are assumed to be constant.

160
A complete description of the equations governing the temporal evolution of fecal pellets and carcasses is provided in the appendices section Text A3 (Eq. A14 and A15). FFGM and GM, and preys are nanophytoplankton, diatoms, microzooplankton, mesozooplankton, small organic particles and large organic particles. A preference of 1 indicates that any prey reached is consumed, a preference of 0 indicates that the prey is never consumed.

Model parameters
Each zooplankton group is characterized by a size range, assuming that sizes within the group are distributed along a spectrum 165 of constant slope -3 in log-log space, according to the hypothesis of Sheldon et al. (1972). These ranges are: 10-200 µm for microzooplankton, 200-2000 µm for mesozooplankton and 2000-20000 µm for macrozooplankton (GM and FFGM).
All parameters in PISCES-FFGM have identical values to those in Aumont et al. (2015). The only exception is the mesozooplankton quadratic mortality rate, whose value has been greatly reduced from its standard value of 3e4 L mol −1 d −1 to 4e3 L mol −1 d −1 since predation by higher trophic levels is now explicitly represented.

170
The parameter values that were introduced in PISCES-FFGM to represent the evolution of GM and FFGM are given in Table 2. Metabolic rates are assumed to vary with size according to the allometric relationship proposed by Hansen et al. (1997).
Therefore, maximum grazing, respiration and flux-feeding rates were calculated from their values for mesozooplankton using a size ratio of 10. The preferences of GM and FFGM for their different prey are detailed in section 2.1.2. Their values are shown in Figure 2. The sinking speed of FFGM carcasses (resp. fecal pellets) is set to 800 m d −1 (resp. 1000 m d −1 ) (Henschke et al.,175 2016). The sinking speeds of GM fecal pellets and carcasses are set rather arbitrarily to 100 m d −1 and 300 m d −1 respectively, within the wide range of values found in the literature (Small et al., 1979;Fowler and Knauer, 1986;Lebrato et al., 2013;Turner, 2015). The quadratic mortality rates have been adjusted by successive simulations evaluated against the observations presented in the next section.

Symbol
Source Table 2. Parameter values used in PISCES-FFGM. The symbols in the "Source" column indicate how the parameter value was determined: ( ) parameters for which we assumed that both GM and FFGM share the same characteristics as mesozooplankton, (•) metabolic rates assumed to vary with size, thus scaled using an allometric scaling convertion of mesozooplankton value based on (Hansen et al., 1997), ( †) parameters tuned to fit PISCES-v2 general biology dynamics, and ( ‡) indicates parameters whose values have been arbitrarily set based on information available in the literature and/or of the authors expertise.

Reference simulation 180
The biogeochemical model is run in an offline mode with dynamical fields identical to those used in Aumont et al. (2015).
These climatological dynamic fields (as well as the input files) can be obtained from the NEMO website (www.nemo-ocean.eu) and were produced using an ORCA2-LIM configuration (Madec, 2008). The spatial resolution is about 2 • by 2 • cos(φ) (where φ is the latitude) with a meridional resolution enhanced at 0.5 • in the equator region. The model has 30 vertical layers with increased vertical thickness from 10 m at the surface to 500 m at 5000 m. PISCES-FFGM was initialized from the quasi-185 steady-state simulation presented in Aumont et al. (2015). The two macrozooplankton groups, their fecal pellets and carcasses were set to a small uniform value of 10 −9 mol C L −1 . The model was then integrated for the equivalent of 600 years, forced with 5-day averaged ocean dynamic fields and with a three-hour integration time step.

Sensitivity experiments
The first experiment, PISCES-GM ("Generic Macrozooplankton"), was designed to investigate the impact of an explicit FFGM 190 representation (with a different grazing parameterization than GM) on the spatial and vertical distribution of POC fluxes. In PISCES-GM, the FFGM ingestion rate (g F F GM m defined in table 1 and used in Eq. A3) was set to 0 which is equivalent to running the model with a single generic macrozooplankton group.
The second experiment, PISCES-HGR ("High Growth Rate"), was designed to investigate the impact of higher clearance rates observed for FFGM than for GM. In PISCES-HGR, the FFGM ingestion rate (g F F GM m defined in table 1 and used in Eq. The third experiment, PISCES-HGM ("High Growth and Mortality rates"), is similar to PISCES-HGR, but tries to compensate the unrealistic high biomasses induced by FFGM high clearance rates in PISCES-HGR. To do so, non-predatory (m c ) quadratic mortality was increased so that FFGM biomass on the upper 100 m is similar to PISCES-FFGM. The quadratic mortality due to predation was not modified because there is no reason to believe that FFGM are subject to a higher predation 200 pressure than GM.
The fourth experiment, PISCES-LOWV ("Low Velocity"), was designed to evaluate the impact of the high sinking speeds of particles from GM and FFGM. In PISCES-LOWV, the sinking speeds of all fecal pellets and carcasses produced by GM and FFGM (w F p X and w Ca X , defined in table 1 and used in Eq. A14 and A15) were assigned the same values as for large particles in PISCES-v2, i.e. 30 m d −1 .

205
The fifth experiment, PISCES-CLG ("Clogging"), was designed to explore the impacts of clogging. Clogging, defined as the saturation of an organism's filtering apparatus with high levels of particulate matter, is a poorly documented mechanism for FFGM but has been observed (Harbison et al., 1986;Perissinotto and Pakhomov, 1997) or suggested (Perissinotto and Pakhomov, 1998;Pakhomov, 2004;Kawaguchi et al., 2004) for some salps species. Unlike other macrozooplankton groups, it has been shown that salps biomass remain relatively low at high chlorophyll concentrations (Heneghan et al., 2020). In 210 PISCES-CLG, the achieved ingestion rate of FFGM (G F F GM , see Eq. A13) is modulated by a clogging function F C (Chl) inspired by the parameterization proposed by Zeldis et al. (1995): In this equation, C th is the clogging threshold, C sh is the clogging sharpness and ERF is the Gauss error function. A low clogging threshold C th of 0.5 mg Chl m −3 is chosen to limit FFGM growth in all moderate and high productivity regions.
Values of the parameters that were changed in the five sensitivity experiments are presented in Table 3. All five sensitivity experiments were initialized with the year 500 output fields from the baseline PISCES-FFGM experiment. They were then run for 100 years. All results presented in this study are average values over the last 20 years of each simulation. PISCES-CLG, PISCES-HGR and PISCES-HGM help to investigate the modeled distribution of GM and FFGM while PISCES-GM 220 and PISCES-LOWV are used for exploring the spatial pattern and depth gradient of particulate organic carbon fluxes.

FFGM biomass estimates
We compiled an exhaustive dataset of in situ pelagic tunicates (i.e., Thaliaceans) concentrations from large scale plankton monitoring programs and previous plankton data compilations to derive monthly field of pelagic tunicates biomass (in mg The Australian Continuous Plankton Recorder (CPR) survey (AusCPR; IMOS (2021)) and the Southern Ocean CPR survey (SO-CPR; (Hosie, 2021)) were excluded because they were found to not quantitatively sample thaliaceans (see appendices, densities in ind m −2 , which we converted to ind m −3 based on the maximum sampling depth of the corresponding net tows.
In KRILLBASE, 5'186 observations of Thaliaceans with missing density values were discarded (35.6% of the original 14'543 of the original 24'316 observations). We examined the composition of the original data sources compiled within JeDI and COPEPOD by assessing the recorded institution codes as well as their corresponding spatio-temporal distributions to evaluate the observations overlapping between these two previous data syntheses. We logically observed a very high overlap between 245 COPEPOD and JeDI as the former data set was the main data contributor to the latter. Therefore, overlapping records were identified based on their sampling metadata, scientific names, concentration values, the recorded institution codes and recorded data sources, and were removed from JeDI. This Most of the records showed a fairly precise taxonomic resolution as 1.6% of the data was species-resolved (mostly S. thompsoni, Soestia zonaria, S. fusiformis and Thalia democratica), 42% genus-resolved (mostly Thalia, Doliolum and Salpa) and 83% family-resolved (mostly Salpidae and Doliolidae). Therefore, we were able to perform taxon-specific conversions 255 from individual concentrations to biomass concentrations (in mg C m −3 ) for each point observation (see Table A1). We used the taxon-specific carbon weights (mg C ind −1 ) summarized by Lucas et al. (2014), which were based on the group-specific length-mass or mass-mass linear and logistic regression equations of Lucas et al. (2011). Not all the observations had a precise counter part in the carbon weights compilation of Lucas et al. (2014) because they were not identified at the species or the genus level (e.g., Class-level, Order-level or Family-level observations). In these cases, we computed the median carbon Hereafter, we will refer to this dataset as "AtlantECO dataset".
Ultimately, monthly Thaliacean biomass fields were computed for validating the monthly FFGM biomass fields of PISCES-270 FFGM. Thaliacea biomass concentrations were averaged per months on a 72x36 grid to obtain the 12 monthly climatological fields of Thaliacea biomass needed for evaluating our model. A low resolution grid (5°x 5°) has been used to counterbalance patchiness of data, as suggested by Lilley et al. (2011). After this final step, the monthly climatological values of Thaliacea biomass concentrations ranged between 0.0 and 454 mg C m −3 (6.53 ± 26.21 mg C m −3 ). Hereafter, we will refer to this climatology as "AtlantECO climatology".

Additional datasets
We also used the monthly fields derived from observations as a standard data set to evaluate some of the other PISCES-FFGM compartments: total macrozooplankton, mesozooplankton, total chlorophyll, nutrients and oxygen.
As with FFGM, for total macrozooplankton observations, a low resolution grid has been used. We use monthly macrozooplankton abundances binned on a 72x36 grid (ind m −3 , vertically integrated between 0 and 100 m) from MARine Ecosystem 280 DATa (MAREDAT) , and then convert abundances to carbon-based concentration to evaluate our modeled distribution of total macrozooplankton biomass (i.e. FFGM and GM). To convert to carbon concentration, an average individual weight of 588 µg was chosen by considering an individual with a mean size of 6.3 mm (the geometric mean of the macrozooplankton size class) and applying the relationship proposed for copepods by Watkins et al. (2011).
Monthly observations fields were binned on a 360x180 grid to validate other modeled distributions. The mesozooplankton

Model evaluation
The model evaluation is based on monthly fields averaged over the last 20 years of the PISCES-FFGM reference simulation.
For each unique observation in the AtlantECO dataset, we sampled the modeled FFGM biomass from the PISCES-FFGM 295 climatology at the corresponding coordinates (latitude,longitude), month, and depth range (minimal depth and maximal depth), so that each observed biomass can be compared to a "model-sampled" biomass. When compared to the AtlantECO climatology, the annual mean FFGM biomass fields and the statistics (Table 4) are calculated from these "model-sampled" biomasses to avoid bias due to sampling.
For other variables, model outputs (N O − 3 , P O 3− 4 , Chl, Mesozooplankton, GM+FFGM) were regridded horizontally and 300 vertically on the same grid as the corresponding observations (see previous section). The macrozooplankton and mesozooplankton fields were integrated vertically on the appropriate vertical range. When compared to observations, model outputs are sampled at exactly the same location and month as the observations. Annually averaged fields as well as statistics (Table 4) are computed from these sampled fields to avoid bias due to sampling. The most striking feature is the reverse distribution of the ratio as compared to the simulated absolute biomass of both GM 315 and FFGM. The ratio exceeds 2 in oligotrophic subtropical gyres while it is minimal in the most productive regions. In eastern boundary upwelling systems, FFGM biomass can be more than two times lower than GM biomass. Vertically, the ratio is on average larger than 1 in the euphotic zone. Below the euphotic zone, it sharply decreases as GM become dominant. In the mesopelagic domain, flux-feeding has been shown to be a very efficient mode of predation (Jackson, 1993 ;Stukel, Ohman, et al., 2019). Since FFGM are not able to practice this feeding mode, they are outcompeted by GM. FFGM:GM ratio is maximum 320 in the lower part of the euphotic zone in the subtropical domain where deep chlorophyll maxima are located. Then, we focus on the evaluation of the new components added in this version of PISCES, i.e. GM and FFGM. In the appendices, we present an evaluation of nitrate, chlorophyll and mesozooplankton (See Text A2 and Fig. A2). For these tracers, note that the performance of PISCES-FFGM is similar to that of PISCES-v2 (Aumont et al., 2015).

325
The annual mean distributions of total macrozooplankton (FFGM and GM) and FFGM only, averaged over the top 300 m of the ocean, are compared to available observations ( Figure 4). A quantitative statistical evaluation of the model performance for these two fields is presented in Table 4. The Spearman correlation coefficient between observed and modeled total macrozooplankton biomasses is 0.26 (p-value < 0.001). Regions of high macrozooplankton biomass are correctly simulated in the northern hemisphere by our model: 94% of the area in which observed concentrations are greater than 0.5 mg C m −3 corre-330 spond to areas in which the simulated concentration is greater than 0.5 mg C m −3 . On the other hand, observations suggest moderate biomass in the Indian Ocean (between 0.05 and 0.5 mg C m −3 ) and low biomass in the Southern Ocean (lower than 0.05 mg C m −3 ). These low and moderate biomasses are not captured by our model, which simulates values greater than 0.5   concentrations are greater than 0.5 mg C m −3 correspond to areas in which modeled concentrations are greater than 0.5 mg C m −3 . However, the maximum observed values are strongly underestimated: the 95th percentile of the modeled values is 2.6 mg C m −3 while it is 32 mg C m −3 in the observations. In the Southern Ocean, the simulated distribution is much more zonally 345 homogeneous than suggested by observations (Fig. 4). Overall, the predicted median biomass of FFGM is similar to that of observations, 0.80 vs. 1.11 mg C m −3 . As with macrozooplankton, but to a lesser extent, the simulated standard deviation is significantly lower than in the observations, 0.96 and 26.9 mg C m −3 respectively. The standard and log10 biases are closer to 0 than those calculated for macrozooplankton (Table 4). We also compared the observed and modeled relationships between FFGM biomass distributions and chlorophyll levels. Here, we present the PISCES-HGR, PISCES-HGM and PISCES-CLG sensitivity experiments and their influence on the FFGM modeled distributions.
A 5-fold increase in the maximum growth rate in the PISCES-HGR experiment leads to a 4-fold and 5-fold increase in the mean and median FFGM concentrations, respectively (see Table 4). While the mean is closer to the observed mean than in the standard experiment, the negative Spearman coefficient shows the unrealistic nature of this simulation and the need 365 to correct mortality accordingly (see Table 4 and Fig. A3). The increase in mortality rate in the PISCES-HGM experiment results in similar mean and median FFGM biomass to the standard PISCES-FFGM experiment (see Table 4) and Fig. A4) but a worse data-model fit (see Table 4 and Fig. A6). Given the large range suggested for the growth rate of FFGM (0.105-1.85 d −1 according to Luo et al. (2022)), these results supports the choice of a conservative approach in our reference experiment (PISCES-FFGM) where the FFGM maximal growth rate is identical to that of GM (i.e. 0.28 d −1 ).

370
The addition of clogging in PISCES-CLG increases the model-data spatial correlation (Spearman's correlation coefficient is 0.34 compared to 0.17 previously, see Table 4  and biases are increased when clogging is added (see Table 4).
None of the sensitivity experiment reproduce the observed spatial variability, which remains much higher than the modeled 380 spatial variability similarly to the standard experiment, and the distribution of observed biomasses is consequently much more spread out than the model (see Fig. A6).

Carbon export from the surface ocean
We first discuss the role of macrozooplankton in shaping the carbon cycle in the upper ocean, focusing on differences between 385 GM and FFGM-related surface processes. Table 5 shows the globally integrated sinking flux of organic carbon particles at 100 m and 1000 m, while Figure 6 focuses on the FFGM-driven carbon fluxes. The total export flux from the upper ocean (at 100 m) is 7.55 Pg C yr −1 (Table 5). This value is relatively similar to previous estimates using different versions of PISCES (Aumont et al., 2015(Aumont et al., , 2017. It is also within the range of published estimates, i.e. 4-12 Pg C yr −1 (e.g., Laws et al., 2000;Dunne et al., 2007;Henson et al., 2011;DeVries and Weber, 2017). Small and large particles produced by phytoplankton, 390 microzooplankton and mesozooplankton account for 91% of this carbon flux. The remaining 9% (0.69 Pg C yr −1 , Table 5) is due to macrozooplankton (FFGM+GM), with one third of this amount coming from carcasses and the remaining from fecal pellets. FFGM are responsible for an export of 0.43 Pg C yr −1 (Table 5), which represents 62% of the total macrozooplankton contribution.
The particularly large contribution from FFGM compared to GM comes from higher production (grazing of 0.94 Pg C yr −1 395 compared to 0.63 Pg C yr −1 for GM, Fig. 6 and S7) while both groups shows similar export efficiency. 45% of the grazed matter is exported at 100 m, the remaining 55% is split between implicit predation by upper trophic levels and loss to dissolved inorganic and organic carbon.

Carbon transfer efficiency to the deep ocean
We then analyze how the representation of the two new macrozooplankton groups influences the fate of particulate organic 400 carbon in the deep ocean. At 1000 m, the total simulated POC flux is 1.97 Pg C yr −1 (Table 5) Table 5. Particulate carbon flux composition at 100 and 1000 m. Units are in Pg C yr −1 . sPOC (resp. bPOC) is for small (resp. large) particulate organic carbon. CaGM (resp. CaF F GM ) is for GM (resp. FFGM) carcasses. F pGM (resp. F pF F GM ) is for GM (resp. FFGM) fecal pellets. from 100 m to 1000 m is 26%. Most of this strong flux reduction is due to the loss of small and large organic particles.
Macrozooplankton-driven export is very effective because it remains almost unchanged from 100 m to 1000 m, 0.69 and 0.67 Pg C yr −1 , respectively (Table 5). Therefore, the contribution of macrozooplankton increases strongly with depth to 34% of the total carbon export at 1000 m (Fig. 7). The respective contribution of particles produced by GM and FFGM (carcasses and 405 fecal pellets) to this flux is almost identical at both depth horizons. At 5000 m, more than 90% of the carbon flux is due to macrozooplankton (Fig. 7).
The PISCES-LOWV sensitivity experiment, in which carcasses and fecal pellets sinking speeds of both macrozooplankton groups are reduced to 30 m d −1 , shows a much greater attenuation of POC fluxes with depth: while the total export of organic carbon at 100 m increases slightly to 7.71 Pg C yr −1 , it is reduced by 20% at 1000 m compared to the standard PISCES-FFGM 410 run (1.56 Pg C yr −1 , see table 5). The macrozooplankton contribution is similar to that found in the standard model at 100 m (8%) but the contribution is reduced to 13% at 1000 m and to 20% at 5000 m (Fig. 7). This confirms that the strong contribution of macrozooplankton to POC fluxes at depth in the standard run is explained by the very high sinking speeds of carcasses and fecal pellets. These high sinking speeds prevent any significant remineralization of these particles as they sink to the seafloor.
The PISCES-GM sensitivity experiment, in which FFGM are not allowed to grow, shows a similar depth gradient of the 415 macrozooplankton contribution (Fig. 7, red curve) compared to the standard run, but a lower contribution at each depth (by 10%). Indeed, the transfer efficiency from 100 to 1000 m differs by only 2% between the two groups in the standard model (97% for FFGM, 95% for GM) so that particles produced at the surface by both groups have a similar fate towards the deep ocean. However, the estimated transfer efficiency is biased as both groups of organisms produce particles below 100 m. Because they can adopt a flux feeding strategy of predation, GM occupy the whole water column whereas FFGM remain confined to 420 the upper ocean (see section 3.1 and Figure 3). As a result, GM also produce particles below 100 m which contribute to the flux at 1000 m and explains the computed higher transfer efficiency. This is confirmed by the PISCES-LOWV experiment: the efficiency of FFGM is reduced to 30% in this simulation while that of GM is only reduced to 40%, even though the carcasses and fecal pellets sinking velocities of both groups are identical. As the remineralization processes are identical in the two runs, we can reasonably assume that the difference comes from the relatively higher productivity below 100 m of GM compared to 425 FFGM.

POC flux spatial patterns
Although the processes underlying the efficient sequestration of the particulate carbon issued from the two groups of macrozooplankton are similar, we investigate how the spatial and temporal patterns of the induced deep POC export differ between GM and FFGM.

430
The relative contribution of FFGM and GM to the POC flux at 1000 m presented in Figure 8 is very contrasted between the two macrozooplankton groups. The POC flux due to FFGM is maximal at about 40% of the total flux in the oligotrophic subtropical gyres. In the productive areas of the low and mid-latitudes, it has intermediate values close to 25%. It is minimal (<15%) at high latitudes, especially along the Antarctic. In contrast, POC fluxes due to GM are maximal in the productive regions of the low and mid-latitudes, especially in boundary upwelling systems where they can exceed 35% of the total flux.

435
These patterns are consistent with the respective spatial distribution of FFGM and GM (ratio shown in figure 3).
We further investigate the importance of GM and FFGM for the spatial patterns of the export of carbon to the deep ocean by contrasting PISCES-FFGM and PISCES-GM experiments (see Section 2.3). Figure 7 shows the relative contribution of macrozooplankton to POC flux as a function of latitude. By comparing the standard model (orange curve) with the experiment without FFGM (PISCES-GM, red curve), we deduce that the explicit representation of FFGM alters strongly the latitudinal 440 distribution of this relative contribution. It is significantly increased at all latitudes. This increase is particularly important in the low latitudes where the contribution goes from less than 20% when FFGM are not allowed to grow (PISCES-GM) to more than 45% in the reference simulation PISCES-FFGM. Furthermore, export due to GM is maximal at about 20°N and S. Compared to GM, the FFGM contribution is relatively constant between these latitudes. This result highlights the strong efficiency of FFGM at exporting organic matter to the deep ocean, in particular in oligotrophic regions with low productivity. The addition 445 of FFGM reduces the contribution of GM at all latitudes, especially at mid and low latitudes in which the contribution losses 15 to 20% (Fig. 7). This reduction results from the competition between FFGM and GM.

Discussion
We added an explicit representation of two macrozooplankton groups in PISCES-FFGM: a generic macrozooplankton group, for which the parameterization is based on an allometric scaling of the mesozooplankton group already existing in PISCES-v2 450 (Aumont et al. (2015), see section 2.2) and which feed mainly on the latter, and an FFGM group that can feed on phytoplankton and microzooplankton. The introduction of FFGM into PISCES, based solely on the representation of their specific diet due to the filter-feeding mode, provided some insights into the potential impacts of FFGM on planktonic communities and carbon cycling at the global scale through trophic effects (e.g. competition with generic macrozooplankton) and efficient carbon export.

Macrozooplankton biomass
After the addition of FFGM in PISCES, our simulation results consistently show that FFGM dominate macrozooplankton in low-productivity regions, but that absolute abundances of FFGM are nonetheless higher in productive areas of the world ocean ( Fig. 3). In a recent study using the COBALTv2 biogeochemical model, Luo et al. (2022) explored the role of pelagic tunicates in the marine ecosystem, with the addition of two new plankton functional groups, i.e. a large salp/doliolid group similar to our 460 FFGM, and a small appendicularian group (Luo et al., 2022). They showed that the FFGM:GM ratio in their model follows a decreasing relationship with chlorophyll, consistently with our modeled FFGM:GM ratio patterns. To better reproduce the relationship between AtlantECO FFGM biomass and chlorophyll from the OC-CCI product, the addition of clogging was needed in our model ( Fig. 5 and section 3.1). Given the paucity of data, it is currently difficult to evaluate these model insights from macrozooplankton databases alone. Heneghan et al. (2020) showed that salps dominate other macrozooplankton groups in 465 low-productivity regions, but, contrary to our model results, these authors also showed that these organisms are more abundant in absolute values in these low-productivity regions than elsewhere in the ocean. Yet, they did not explore the processes that could drive this distribution. As evidenced by our PISCES-CLG experiment, clogging may be a potential explanatory mechanism but the evidence for this process is weak. Future studies are needed to determine the processes involved in limiting FFGM biomass at high chlorophyll concentrations.

Export of organic carbon
Our modeled FFGM have a weak impact on phytoplankton and microzooplankton biomasses, due to the low predation pressure they exert on these low-trophic levels (grazing flux of 1 Pg C yr −1 , which represents less than 3% of primary productivity).
Nevertheless, due to the high sinking speed of FFGM-derived fecal pellets and carcasses, FFGM substantially increase the carbon export ratio and transfer efficiency. We compiled results from distinct studies on global biogeochemical impacts of 475 FFGM in table 6 to support our results.
The overall PISCES-FFGM modeled production of POC by FFGM in the upper 100 m is 0.42 Pg C yr −1 (Table 6). This value falls within the range of data-driven estimates (Table 6). It is an order of magnitude above the value of 0.03 Pg C yr −1 from Lebrato et al. (2019), presented as a lower bound estimate due to their conservative assumption of equivalence between GZ annual production and total GZ biomass. On the other hand, our simulated FFGM POC production within the top 100 480 m is 10 times lower than the estimate of 3.9 Pg C by Luo et al. (2020). In this study, FFGM production was forced offline by modeled phytoplankton and zooplankton climatologies, so that FFGM predation had no feedback on their prey biomass. Luo et al. (2020)'s production estimate can be seen as an upper estimate as GZ-induced predation pressure would affect the biomass of other trophic levels in a fully-coupled model, thus affecting the gelatinous biomass itself and the induced carbon fluxes. Indeed the higher FFGM POC production is mostly due to a higher FFGM grazing in their study (6.6 Pg C yr −1 485 compared to our modeled value of 1 Pg C yr −1 , Table 6). Finally, our modeled FFGM impacts on upper ocean POC are similar to those by Luo et al. (2022) based on COBALT-GZ: the simulated production of detritus by FFGM in the first 100 m in our model is twice lower than in Luo et al. (2022) and the effective export of these detritus at 100 m is 30% lower ( Table 6). The smaller difference in export than in production lies in the use of a 10-times lower particle sinking speed and a 20-times higher remineralization rate in COBALT-GZ (Stock et al., 2014) compared to PISCES-FFGM, resulting in a lower production export   is for export to. Teff is for transfer efficiency. Tot MAC is for total macrozooplankton (GM + FFGM). * Lebrato et al. (2019)  detritus in the upper 100 m than large tunicates, which supports our choice to represent only FFGM (i.e. macrozooplankton) and not filter-feeding mesozooplankton in our biogeochemical model.
The impact of an explicit representation of FFGM on POC export is negligible in both models when compared to a version without FFGM (+/-2%, Table 6). But the contribution of total macrozooplankton to POC fluxes increases significantly with 495 FFGM in both models (GZ-COBALT: +41%, PISCES-FFGM: +55%, Table 6) and this despite the simulated decrease in export by GM (-11% in GZ-COBALT, -19% in PISCES-FFGM, Table 6), so that the contribution of FFGM only to POC export at 100 m in both models is more than 5% (Table 6). Thus, we can reasonably state that the representation of FFGM in a biogeochemical model redistributes the carbon particles between the different compartments over the top 100 m (more of very large particles from macrozooplankton, less of small particles from smaller organisms) without significantly altering the 500 total amount. This change in particles composition is key to the major role that FFGM play in the export of carbon to the deep ocean.

Deep carbon fluxes
FFGM have a modest impact on subsurface export (less than 10 % of the global POC export at 100 m depth), but this impact is highly increasing with depth, reaching much higher values at the seafloor (>40%) and suggesting that FFGM play a key role 505 in carbon sequestration in the deep ocean. We also demonstrated that surface FFGM productivity and the transfer efficiency of FFGM-driven POC are key processes that strongly affect the magnitude and distribution of deep POC export.
The FFGM-driven export of POC at 1000 m (resp. seafloor) of 0.42 (resp. 0.39) Pg C yr −1 falls between the low value of 0.02 (resp. 0.01) Pg C yr −1 ) proposed by Lebrato et al. (2019) and the much larger estimate of 1.4 (resp. 0.86) Pg C yr −1 given by Luo et al. (2020) (Table 6). The quite large differences between these estimates are mainly explained by the evaluation of In addition to surface productivity, the efficiency of POC transfer is critical to the absolute value of POC export at depth.
The sinking velocity of particles is a key factor that strongly controls this efficiency. In the studies of Lebrato et al. (2019) and 515 Luo et al. (2020), in which the sinking velocities are greater than 650 m d −1 , the transfer efficiency is about 50% (Table 6). It is reduced to 25% when the FFGM fecal pellets (which account for 80% of FFGM detritus in their study) velocity is reduced to 100 m d −1 in Luo et al. (2020). The same finding was obtained when reducing the velocity from 800-1000 m d −1 to 30 m d −1 in our experiment PISCES-LOWV, where the transfer efficiency from 100 to 1000 m decreases from 97% to 30%. However, due to the use of a low remineralization rate, our simulated transfer efficiency from 100 to 1000 m is very high compared to 520 Luo et al. (2020) for similar carcasses and fecal pellets sinking speeds (Table 6). Still, our transfer efficiency in PISCES-FFGM fits the vertical profiles of depth attenuation of jelly-driven organic matter export proposed by Lebrato et al. (2011) for high sinking velocities and low remineralization rates.
Last but not least, PISCES-FFGM seems to capture the intensity and part of the variability of the intense carbon export events described by Henschke et al. (2016) linked to short time proliferation events of FFGM: they estimated the export potential at 525 species, the minimum from 0.6 to 1171 mg C m −2 and the maximum from 656 to 77 143 mg C m −2 . We compare these results to the annual maxima of the FFGM carbon export simulated at each grid point by our model (Table 6). The values obtained range from 0.34 to 1580 mg C m −2 with a spatial mean of 141 mg C m −2 , which is consistent with the species-range of mean, min and max in their study (Table 6). This also supports our choice of a very low remineralization rate and high sinking rates.

530
The latter is confirmed with the PISCES-LOWV experiments in which modeled export maxima fall below the min, mean and max ranges of Henschke et al. (2016).

Data-based climatology
To evaluate the modeled FFGM biomasses, we compiled data from different sources (section 2.4) to produce a gridded climatology of large pelagic tunicates. Our AtlantECO dataset is based on similar observations as the previously compiled dataset 535 (Luo et al., 2020(Luo et al., , 2022), but we used a different approach to convert abundances to biomasses by taking into account the taxonomic information available on the samples, even when the species is not indicated.
Our model predicts a median biomass of FFGM similar to our dataset (0.80 vs. 1.11 mg C m −3 ), and reproduces 91% of the areas where biomass is high (>0.5) ( Table 4). The introduction of a clogging mechanism, which would represent a saturation of the salp filtering apparatus for high prey concentrations, improves the representation of low biomass areas (section 2).

540
In PISCES-CLG, a sensitivity experiment in which the clearance rate is decreased for chlorophyll concentrations above 0.5 mg Chl m −3 , the Spearman correlation coefficient is doubled when comparing simulated and observed FFGM concentrations.
Note however that this clogging mechanism and its impact on pelagic tunicates growth are largely under-documented, and rely on a few 30-yr old publications (Harbison et al., 1986;Fortier et al., 1994).
However, our modeled variability of the spatial distribution of FFGM was 25 times lower than the observed variability 545 (Table 4). This large variability in observations has already been described in previous compilations of pelagic tunicates observations (Luo et al., 2020(Luo et al., , 2022. Numerous aspects may contribute to the high variability of observations compared to models: scarcity of the observations, design of the sampling strategy (Hjøllo et al., 2021), biases in the sampling and enumeration methods (Frank, 1988;Mack et al., 2012), use of species-and location-dependent conversion factors (Arhonditsis and Brett, 2004), differing definitions of the compared groups or communities and the scale of investigation (local measure-550 ments are compared to average 5x5°estimates). Indeed, zooplankton patchiness increases with organism size . Physical (mesoscale and submesoscale processes) and biological (diel vertical migrations, predator avoidance, food patches, mate search) processes combine to drive zooplankton patchiness (Folt and Burns, 1999). Although the introduction of a macrozooplankton compartment (namely cnidarian jellyfish) has been shown to increase patchiness in a recent modeling study (Wright et al., 2021), the spatial resolution (≈2 degrees) of our model setup, and the lack of key biological processes 555 (e.g., complex life cycle) in our model likely preclude the representation of such patchiness. Another source of uncertainty lies in the use of a taxon-specific carbon conversion factor to convert thaliacean abundance data to biomass data. While this approach is appropriate for many protists, thaliacean biomasses estimates based on this method are highly uncertain because these organisms can vary in length by more than one order of magnitude (Iguchi and Ikeda, 2004). In particular, most of the time, when a net returns hundreds of salps, these salps are relatively young blastozooids (i.e., on the small end of the size 560 range). Thus estimating biomass from abundance may lead to an overestimation of the true biomass variability. This supports the need to move towards a systematic reporting of biomass (or at least biovolumes) during zooplankton surveys.
Also, the data temporal resolution is insufficient to analyse seasonal patterns: only 7% of the grid points in the AtlantECO climatology are derived from data covering at least 6 distinct months. Yet, our standard PISCES-FFGM simulation shows an approximate one-month lead in the seasonal biomass peak of FFGM compared to GM, this lag being consistent at the global 565 scale to that of the food of the two groups ( Figure A8). This suggests that the filter-feeding mode of FFGM may have an impact on the temporal dynamics of the FFGM-driven POC flux. However, it is difficult to give a high confidence level to this statement because the spatial distributions between the lags of the organisms and their food are very patchy and the temporal variability of the prey does not correspond to that of the corresponding groups when focusing on specific regions ( Figure A8).
This claim supports the need to improve the temporal monitoring of FFGM populations in order to understand their seasonality 570 and thus characterize the seasonal variations of FFGM impacts on carbon fluxes.

Boom-and-burst dynamics
Pelagic tunicates exhibit pullulation-extinction population dynamics, i.e. the alternation between rapid growth phases and massive mortality events. As a consequence, patchiness is particularly strong for gelatinous zooplankton (Graham et al., 2001;575 Purcell, 2009;Lilley et al., 2011;Lucas et al., 2014). However, this dynamics is clearly not simulated by PISCES-FFGM.
This result was expected as biogeochemical models are known to struggle to reproduce the observed spatial variability in the abundance of different groups of meso-and macro-zooplanktonic organisms (Wright et al., 2021). From a biogeochemical perspective, the impacts of FFGM on ecosystem structure and carbon export are therefore "smoothed" in time and space when simulated by PISCES-FFGM. Still, the results obtained provide a first assessment of the annual impacts of FFGM at the global 580 scale and in large biogeochemical regions (e.g. low productive oligotrophic gyres vs. highly productive upwelling regions).
However, the currently modeled FFGM ability to consume prey over a wide size range is not the only factor likely to trigger boom-and-bust dynamics. FFGM high clearance rates and complex life cycles with an asexual reproductive phase, currently not represented in the standard model, are also likely to play a role in such dynamics. In the PISCES-HGR sensitivity experiment, increasing growth rates of FFGM without adequate modifications of FFGM mortality rates caused the generic 585 macrozooplankton population to collapse because they were outcompeted by FFGM everywhere except in the mesopelagic and deep ocean. As expected, and similarly to Luo et al. (2022), the modification of the quadratic mortality in the PISCES-HGM sensitivity experiment neither improved the fit with the observations, nor triggered any boom-and-burst dynamics. To further investigate the effect of high growth rates and clearance rates of FFGM, a better understanding of the physiological and environmental drivers of the FFGM mortality processes triggering the end of their swarms seems essential, as their causes are 590 multiple and too poorly documented to be currently modeled .
Also, life cycles are currently not represented in the model though it could significantly affect the temporal dynamics of a biogeochemical-model (Clerc et al., 2021). Most FFGM have a complex life-cycle, with an alternation between a sexual and asexual phase that could be a major driver of their population dynamics (Henschke et al., 2016). A single-species observation based study on Thalia democratica in South-East Australia suggested that life history characteristics such as asexual reproduc-595 tion and growth are associated with inter-annual variations in abundance and thus may be major factors determining population dynamics, in particular the magnitude of swarms (Henschke et al., 2014). Inclusion of such life cycle traits in a single-species model of Salpa thompsoni in the Southern Ocean helped understand the seasonal and interannual variability of salp abundance (Henschke et al., 2018). These studies are focused on one species and one region, and the inclusion of their life cycle in a global model in which FFGM constitute a single compartment would require a multispecies large scale evaluation of the FFGM life 600 cycle role in the temporal dynamics of the swarming process.

Carcasses and fecal pellets transfer efficiency
One of the greatest sources of uncertainty about the export of carbon from FFGM to the deep ocean is the transfer efficiency (see Table 6), which depends primarily on remineralization rates and sinking speeds. This raises questions about the processes that could affect the fate of carcasses and fecal pellets (CAFP) as they sink. At a given temperature, our simple FFGM representation 605 includes constant remineralization of CAFP and consumption through filter feeding by GM (Eq. A14 and A15). The induced losses are very low compared to FFGM's CAFP production rates (<5%). However, predation by scavengers could significantly affect CAFP during their fall (Dunlop et al., 2018;Scheer et al., 2022). Benthic consumption by scavengers is well documented for jellyfish carcasses (Sweetman et al., 2014;Henschke et al., 2013), but their fate in the vertical column is largely unknown.
Also, parasitism by hyperiid amphipods is likely to affect FFGM carcasses production and degradation, and thus affects deep 610 carbon export by FFGM (Lavaniegos and Ohman, 1998;Phleger et al., 2000;Hereu et al., 2020). Lastly, most measured sinking speed values are based on small (a few meters) sinking column experimental setup and thus do not account for any degradation process (Lebrato et al., 2013). Thus, by combining particularly high velocities with a partial representation of the degradation processes, we mechanistically obtain a particularly high transfer efficiency of FFGM particles. Our estimate of the impact of FFGM on the deep carbon cycle should therefore be interpreted as an upper bound, and a better understanding of FFGM 615 carcasses and fecal pellets fate is needed to properly estimate their impacts on the deep ocean.

Conclusions
We explicitly represented large pelagic tunicates in the global marine biogeochemistry model PISCES and evaluated the simulated distribution of FFGM by compiling available observations into a FFGM biomass climatology using a taxon-resolving biomass-abundance conversion. Representation of FFGM in a marine biogeochemical model has a small impact on total de-620 tritus production in the first 100 m, with 6% of this production due to FFGM. Due to their high sinking speeds, almost all of the organic matter produced by FFGM is transferred to the deep ocean. Therefore, FFGM carcasses and fecal pellets dominate the export of organic matter in the deep ocean (e.g. 70% at 5000 m). The spatial distribution of FFGM-driven export differs from that of the other macrozooplankton group, GM, which also contributes significantly to export at depth (25% at 5000 m). Indeed, due to their filter-feeding mode of predation, access to preys of variable size allows FFGM to better exploit low 625 productivity environments than GM, especially in subtropical oligotrophic gyres, where FFGM are twice as abundant as GM and thus contribute 5 times more to POC export at 1000 m.
A more detailed inclusion of the processes involved in the boom-and-burst dynamics of FFGM (e.g. life cycle, clogging, high clearance rates) will be necessary to better understand the spatial and temporal variability of their impacts on carbon export and ecosystem structure. Still, a promising perspective would be to run our PISCES-FFGM model forced by climate 630 projections. Such a simulation would allow analysis of annual global and large scale regional trends in the impact of FFGM on marine biogeochemistry. In particular, as climate change could favor small phytoplankton (Peter and Sommer, 2013), we could expect an amplification of the spatial pattern we currently described, with FFGM even more favored in low productive regions.
Code and data availability. This section needs to be completed. All raw and gridded data sets will be made publicly available in open access 635 within the framework of the European H2020 project AtlantECO (grant agreement no 862923). Preliminary DOIs can be made available to the reviewers upon request. All model outputs necessary to reproduce the results in this manuscript will be made publicly available. 3. Table A1 Text A1.
When including Aus-CPR and SO-CPR data, the resulting point biomass measurements ranged between 0.0 mg C m −3 and 19'451 mg C m −3 , with and average of 0.63 ± 48 mg C m −3 . However, this range is largely zero-inflated (94.6% of the 645 observations corresponded to a biomass of 0.0 mg C m −3 ) due to the high relative contribution of both CPR surveys whose data only comprised 1.1% of non null values. Such strong zero inflation can be attributed to sampling artifacts due to the specificities of the CPR and thus very likely do not reflect reals absences (Richardson et al., 2006). Indeed, the CPR continuously collects plankton at standard depth of 7 m and at a speed of nearly 0.2 m s −1 , as seawater flows in through a square aperture of 1.61 cm 2 , which is too narrow to adequately sample large gelatinous macrozooplanton such as salps and doliolids, especially in 650 the Southern Ocean (Pinkerton et al., 2020). Consequently, we decided to remove the observations from the AusCPR and the SO-CPR from our final validation data set.

Nutrients
Map a. (resp. b.) in Fig. A2 presents the observed (resp. simulated) surface concentrations of nitrates. The model performs 655 particularly well for surface nitrates, with absolute values and simulated spatial patterns very consistent with observations (r=0.83). The model performance is very similar for phosphates (r=0.83).

Chlorophyll
The modeled annual chlorophyll distribution is compared to OC-CCI satellite observations in Fig. A2 c. and d. The correspondence between the observed and simulated surface chlorophyll is rather satisfactory (r= 0.59). The average value is similar  of about 2 to 2.5 by the algorithms deducing chlorophyll concentrations from reflectance as discussed in Aumont et al. (2015).

Mesozooplankton
Mesozooplankton annual distribution on the top 300 m is compared to the MAREDAT product in Fig. A2 e. and f. The model performs quite well (r=0.45) and fits the observed spatial patterns, and the distribution of high vs. low concentration regions.
However, it tends to overestimate the low concentrations and underestimate the high concentrations. Indeed, mesozooplankton variability is slightly reduced in the model (standard deviation of 0.34 vs 0.59 mmol C m −3 in the observation).

Macrozooplankton dynamics:
G X , the ingested matter, is depending on food availability to X. We distinguish two predation behaviours: concentrationdependent grazing and flux feeding.
Concentration-dependent grazing is based on a Michaelis-Menten parameterization with no switching and a threshold (Gen-675 tleman et al. , 2003). The equation describing the grazing rate of X on prey I, g X (I), is derived as: where F X is the available food to X, g X m is the maximal grazing by X rate, F X thresh is the feeding threshold for X, I X thresh is 680 the group I threshold for X, K X G is the half saturation constant for grazing by X, p X I is the X preference for group I. Flux-feeding accounts for particles traps deployed by some zooplankton species (Jackson, 1993). It is derived as a particles flux depending term, an thus depends on the product of the concentration by the sinking speed: where ff H (I) is the flux-feeding rate of prey X on particle I, ff H (I) is the maximal flux-feeding rate of prey X on particle 685 I, w I is the vertical sinking velocity of I particles.
For GM: with E ff GM the proportion of filter-feeders, G maxff GM the potential ingestion by flux feeding,G ff GM the actual ingestion by flux feeding , G g GM the ingestion by concentration dependent grazing and p X Y the X preference for group Y

695
For FFGM: For the PISCES-CLG experiment (with FFGM clogging) run, the ingested matter by FFGM G CLG F F GM is: where F C (Chl) is the clogging function presented in Eq. 2 of the paper.

Carcasses dynamics:
Carcasses production by organisms X (=F F GM or =GM ) comes from non predatory quadratic and linear X mortalities.
Loss terms include a temperature dependent term representing remineralization by saprophagous organisms and flux-feeding by GM. Flux feeding includes two terms : the ingested food by GM which is temperature dependent and the non ingested 705 matter fractionated by flux feeding process (Dilling and Alldredge, 2000), which is assumed to be equal to the ingested portion except the temperature dependency.
Where α is the remineralization rate.

Fecal pellets dynamics:
Fecal pellets production by organisms X (=F F GM or =GM ) comes from non assimilated food. Loss terms, similarly to 715 carcasses, include a temperature dependent remineralization term and a flux-feeding by GM term.       of the diagram represents the inflows and outflows of GMs integrated globally over the first 100 meters. The inflow is the grazing on the different prey. The arrow going from GM to GM corresponds to the flux related to growth due to assimilated food. The outflows are : i) the remineralization/non-assimilation processes that go into the dissolved organic carbon (DOC) and dissolved inorganic carbon (DIC) ii) the quadratic and linear mortality terms (directly remineralised in PISCES-FFGM because of the lack of explicit representation of upper level predators) and iii) the production of particular organic carbon (POC) via carcasses and fecal pellets. The lower part of the diagram corresponds to the export of POC linked to the fall of carcasses and fecal pellets of GM. The values in blue correspond to the global annual GM-driven POC flux through the corresponding depth, the values in parenthesis representing the total POC flux (related to FFGM, GM, bPOC and sPOC). Figure A8. Spatial distribution of the annual period of maximum macrozooplankton biomasses and maximum food availability A filter was applied to keep only the areas at more than 20°latitude from the equator and in which the amplitude of annual biomass variation is higher than 20%. The amplitude is calculated as (2 × (max − mix)/(min + max)) with min the minimum annual biomass and max the maximum annual biomass. a. Map of months with maximal FFGM biomasses b. Map of lag (in months) between months of maximal FFGM biomasses and months of maximal FFGM biomasses c. Map of months with maximal FFGM food availability (calculated as the sum of prey weighted by FFGM preferences for each prey) d. Map of lag (in months) between months with maximal FFGM food availability and months with maximal GM food availability.