Bridging the gaps between particulate backscattering measurements and modeled particulate organic carbon in the ocean

Oceanic particulate organic carbon (POC) is a relatively small (~4 Pg C) but dynamic component of the global carbon cycle with fast mean turnover rates compared to other oceanic, continental and atmospheric carbon stocks. Biogeochemical models historically focused on reproducing the sinking flux of POC driven by large fast-sinking particles 10 (bPOC). However, suspended and slow-sinking particles (sPOC) typically represent 80–90% of the POC stock, and can make important seasonal contributions to vertical fluxes through the mesopelagic layer (200–1000 m). Recent developments in the parameterization of POC reactivity in the PISCES model (PISCESv2_RC) have greatly improved its ability to capture sPOC dynamics. Here we evaluated this model by matching 3D and 1D simulations with BGC-Argo and satellite observations in globally representative ocean biomes, building on a refined scheme for converting particulate backscattering 15 profiles measured by BGC-Argo floats to POC. We show that PISCES captures the major features of sPOC and bPOC as seen by BGC-Argo floats across a range of spatiotemporal scales, from highly resolved profile time series to biomeaggregated climatological profiles. Our results also illustrate how the comparison between the model and observations is hampered by (1) the uncertainties in empirical POC estimation, (2) the imperfect correspondence between modelled and observed variables, and (3) the bias between modelled and observed physics. Despite these limitations, we identified 20 characteristic patterns of model-observations misfits in the mesopelagic layer of subpolar and subtropical gyres. These misfits likely result from both suboptimal model parameters and model equations themselves, pointing to the need to improve the model representation of processes with a critical influence on POC dynamics, such as sinking, remineralization, (dis)aggregation and zooplankton activity. Beyond model evaluation results, our analysis identified inconsistencies between current estimates of POC from satellite and BGC-Argo data, as well as POC partitioning into phytoplankton, heterotrophs 25 and detritus deduced from in situ bio-optical data. Our approach can help constrain POC stocks, and ultimately budgets, in the epipelagic and mesopelagic ocean. https://doi.org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction 30
The biological carbon pump (BCP) is the ensemble of processes that transfer the organic matter produced by plankton in the sunlit ocean surface to deeper layers (Volk and Hoffert, 1985). This vertical flux plays a central role in the Earth's climate, as it influences the oceans' capacity to absorb and ultimately store atmospheric CO 2 over centennial or millennial time scales (Kwon et al., 2009). The BCP is also central to biogeochemical functioning of the ocean, as it determines the quality and quantity of organic matter available to the ocean interior Hernández-León et al., 2020) and the 35 interior . BGC-Argo floats provide vertical profiles of temperature, salinity, bio-optics and chemical variables between 0-1000 m every 1 to 10 days in near-real time, and are thus especially well suited to study particles in the 95 mesopelagic layer, where the strongest POC gradient occurs. The rapidly growing fleet of BGC-Argo floats equipped with bio-optical sensors enables comparison between models and observations at global scales with enhanced spatiotemporal resolution. Unfortunately, BGC-Argo floats measure only a bio-optical proxy of POC, the particulate backscattering coefficient (usually at 700 nm, b bp700 ) and empirical conversion factors are needed to estimate POC from b bp700 (Cetinić et al., 2012;Stramski et al., 2008). 100 In this study we compare sPOC and bPOC concentrations estimated from BGC-Argo floats to their PISCES-simulated counterparts, as well as satellite-retrieved surface POC. The comparison is enabled by a novel empirical algorithm to convert b bp700 to POC. Observations and simulations are matched in 3D (biome-wide climatological scale) and 1D (at defined locations over an annual cycle). These complementary strategies allow us to evaluate the skill of PISCES at simulating POC 105 stocks and fractions, and indirectly the drivers of POC export and mesopelagic POC budgets, in globally representative biomes.

Definition of vertical and horizontal domains
Studies of the BCP usually decompose the ocean into vertical domains: a surface layer where autotrophic activities 110 dominate, and one or several ocean interior layers where heterotrophic processes dominate. Functional definitions based on light penetration, peak export production, vertical mixing or long-term carbon sequestration are usually the most appropriate ones for process studies (Buesseler and Boyd, 2009;Buesseler et al., 2020;Guidi et al., 2015;Palevsky and Doney, 2018).
Because this paper is mainly descriptive and combines observations and simulations, we will refer to layers defined by fixed depths: epipelagic (0-200 m), mesopelagic (200-1000 m) and bathypelagic (1000-4000 m). 115 Over the horizontal dimensions, our comparisons between observations and model results rely on the ocean biomes defined by Fay and McKinley (2014). These authors subdivided each ocean basin (Atlantic, Pacific, Indian, and Southern Ocean) into different biomes on the basis of observed variables, namely: sea surface temperature, spring/summer chlorophyll a concentrations (Chla), ice fraction, and maximum mixed layer depth (MLD), all on a 1°×1° grid. This division resulted in 17 120 regions ascribed to one of the following five biomes: the ice biome, the subpolar seasonally stratified biome, the subtropical seasonally stratified biome, the subtropical permanently stratified biome, and the equatorial biome. The analyses reported herein focus on the following four biomes (Fig. 1): the seasonally stratified North Atlantic subpolar gyre (NASPG); the permanently stratified Atlantic and Pacific subtropical gyres, which were grouped together (STG); the seasonally stratified Southern Ocean (Subantarctic); and the Mediterranean Sea, which was added here owing to the abundance of BGC-Argo 125 https://doi. org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c Author(s) 2021. CC BY 4.0 License. data, and represents a seasonally stratified subtropical biome. Fay and McKinley's definition allows biome boundaries to change from one year to another. Here we analysed only data from the core of each biome, defined as the grid cells that never changed classification during the 1998-2010 satellite observation periods.

BGC-Argo observations
The global dataset acquired by the array of BGC-Argo floats was downloaded on 14 January 2020 from the Global Data 130 Assembly Center hosted by Ifremer (ftp://ftp.ifremer.fr/ifremer/argo/dac/) (Argo, 2000). We selected floats equipped with at least a Chla fluorometer and a sensor for optical backscattering at 700 nm (b bp700 ), in addition to the conductivitytemperature-depth (CTD) probe. The downloaded measurements had undergone the standard processing, which includes the application of calibration equations to raw sensor output and the performance of near-real-time quality control to both CTD (Wong et al., 2021) and Chla measurements (Schmechtig et al., 2018). Since no specific quality control procedure has been 135 established yet for b bp700 profiles, the latter were only flagged according to the general criteria (Schmechtig et al., 2016).
Thus, we used all b bp700 measurements with quality control flag ≤ 3 (equivalent results were obtained with flag ≤ 2). Two different processing pipelines were applied to different subsets of the BGC-Argo data, as described below.

Global gridded climatologies (3D approach)
The global dataset acquired between 2010 and 2019 was used to produce global gridded monthly and seasonal climatologies 140 for b bp700 and Chla. The measurements were binned onto the ORCA2_L31 grid used for NEMO-PISCES simulations (see 2.4.1), which has an horizontal resolution of about 2° that increases to 0.5° in the meridional direction in the equatorial domain, and 30 oceanic vertical levels between the surface and the ocean bottom. The thickness of the vertical bins increases progressively from 10 m at the surface to 339 m in the 22 nd bin (870-1209 m), the deepest one containing BGC-Argo data.
In each grid element, the average, median, range and data counts were computed. Profiles from the CSIRO and INCOIS data 145 assembly centers were not used because, at the time of download, they had not taken into consideration the new calibration files provided by the manufacturer. A total of 72460 profiles were used to calculate the global gridded climatologies. the σ θ at 5 m after applying a 5-point running mean to the top 10 m of the profile. Eleven additional MLD criteria were also calculated to assess the robustness of the approach (Fig. S1). Following Briggs et al. (2011Briggs et al. ( , 2020, each b bp700 vertical profile was smoothed with sequential 11-point running-minimum and running-maximum filters to separate the baseline from the spikes. The baseline signal corresponds to the bulk population of small particles, whose diameter is smaller than 100 µm and 160 mostly between 0.5 and 30 µm (Dall'Olmo et al., 2009;Organelli et al., 2018). Each spike reflects the passage of a particle (organism or aggregate) larger than about 100 µm in front of the sensor window. Unlike Briggs et al. (2020), we did not subtract from the baseline profile the 850-900 m signal, which in that study was attributed to a background of small refractory particles with constant concentration. The baseline and spike signals were converted to sPOC and bPOC, respectively, as described in the next section. All measurements were subsequently averaged into 18 vertical bins of 165 progressively increasing thickness, such that the deepest bins contained at least 10 measurements. Finally, each profile was interpolated onto the L75 vertical grid commonly used in NEMO simulations. This grid has 46 bins between the surface (0-1 m) and the deepest layer considered here (901-996 m). The profile time series was temporally binned into 5-day periods.
For the comparison to PISCES 1D simulations, BGC-Argo time series were cut into one-year periods (shifted by 6 months in 170 the Southern hemisphere), which we will call coherent annual time series (CATS) hereafter. The CATS fulfilled the following conditions: (1) sampling dates spanned at least between days of year 25 and 340; (2) the float remained in the same region and did not cross major oceanic fronts according to the vertical-temporal evolution of temperature, salinity, σ θ and spiciness; (3) bottom depth was > 1000 m for all profiles (bathymetry obtained from the 15 arc-second GEBCO 2019 product; https://www.gebco.net/data_and_products/gridded_bathymetry_data/gebco_2019/gebco_2019_info.html); and (iv) 175 the Chla and b bp700 sensors were stable according to both the vertical profiles and the continuous measurements acquired during drift at 1000 m between profiles. A total of 50 CATS from 28 different floats were selected, with 10-16 CATS in each biome, 32 (18) in the Northern (Southern) hemisphere (Table S1).

Conversion of b bp700 to POC
To convert the profiles of the backscattering coefficient at 700 nm (b bp700 ) to POC we developed an empirical algorithm, 180 whose behaviour is summarized in Fig. 2, building on previous studies Evers-King et al., 2017). The variability of POC/b bp700 enclosed in this algorithm results from processes that alter the size distribution, dominant shapes, and chemical composition of the particle assemblage, and ultimately its bulk optical properties, as discussed in section 4.2.
Further details are provided in the Appendix. The algorithm estimates the POC/b bp700 ratio along the vertical profile between 0 and 1000 m, and proceeds in two steps. First, the POC/b bp700 ratio is calculated by prescribing a POC/b bp700 ratio in the 185 surface layer (z surf,biome ) and an exponential decrease with depth in the underlying water column, which converges asymptotically towards a constant deep value (c): POC/b bp,700 (z) = c + a biome *exp(-b*(z-z surf,biome )*0.001), z > z surf,biome (1)

190
The a biome coefficient is biome-specific, whereas the asymptote at depth is fixed at c = 1000 mmol C m -3 m. The z surf,biome corresponds to the 5% quantile of the climatological MLD in summer in a given biome. The POC/b bp,700 ratios at z surf,biome , corresponding to a biome + c in Eq. 1, are taken from the literature and range between 2600 and 4900 mmol C m -3 m ( Fig. 2; The exponential decrease prescribed by Eq. 1 is similar to that proposed by Bol et al. 2018, except for the inclusion of the 200 constant term c that prevents the ratio from becoming 0 at depth. The slope of the exponential decrease (b = -6.57) is constant in all biomes and based on our fit to the Cetinić et al. (2012) dataset, using the same depth bins as Bol et al. (2018), but additionally forcing the curve towards c at 1000 m.
The uncertainty of regional b bp700 -POC conversion factors in the epipelagic is typically <10% according to the standard error 205 of the POC vs. b bp700 linear regression slopes (Table A1). The few available measurements in the mesopelagic suggest a POC/b bp700 uncertainty lower than a factor of two. Through this study, we will assume that model/observation ratios larger/smaller than 2/0.5 can safely be regarded as model over/underestimates, which possibly is a conservative criterion for the epipelagic layer.

210
The conversion of b bp700 to POC was done using different MLD data for the global climatologies and the CATS. For the global climatologies we used the Monthly Isopycnal & Mixed-layer Ocean Climatology (MIMOC) of Schmidtko et al.
Although MIMOC is based on an algorithm that evaluates several MLD criteria, it has been shown to be globally consistent with the MLD based on a 0.03 kg m -3 σ θ threshold (Holte and Talley, 2009;Sallée et al., 2021). For the float time series, we 215 used the MLD defined by a 0.03 kg m -3 threshold computed for each individual profile.

Ocean color satellite data
Satellite observations for the 1997-2019 period were downloaded from GlobColour (https://www.globcolour.info), a merged multisensor dataset, on 2 March 2020. Monthly sea-surface POC fields based on the algorithm of Stramski et al. (2008) were used to compute monthly climatologies that were subsequently reprojected onto the ORCA2 grid.

PISCES simulations and matching with observations
Simulations were run using the ocean biogeochemistry model PISCESv2 (Aumont et al., 2015) with the RC parameterization for detrital POC (Aumont et al., 2017). The configuration of PISCES used here has 24 tracers: two classes of phytoplankton ("nanophytoplankton" and diatoms), detrital particles (small and big), and zooplankton (micro-and mesozooplankton), plus 18 additional tracers that comprise dissolved inorganic macronutrients and iron, inorganic carbon chemistry variables, 225 dissolved organic carbon (DOC), and different particulate stocks of iron and silica. Phytoplankton growth depends on light, inorganic nitrogen, phosphorus and iron, with an additional silicate requirement for diatoms. Microzooplankton and mesozooplankton consume the two classes of phytoplankton and detrital particles with different preferences, and mesozooplankton also predate on microzooplankton. Detrital particles are produced through the mortality of phytoplankton and zooplankton (which are routed to small and big particles in different proportions), zooplankton sloppy feeding, and DOC 230 coagulation. Production of large detritus also results from enhanced diatom mortality upon bloom collapse, aggregation of small detritus, and zooplankton mortality and fecal pellet production (the latter two derived from a closure term that accounts for unresolved higher trophic levels). Small and big detrital particles are nominally smaller/larger than 100 µm and sink, respectively, at 2 and 50 m d -1 . Both small and big detritus are remineralized by implicit bacterial activity, and consumed by flux-feeding mesozooplankton. Remineralization follows first-order kinetics with an initial specific rate "k" of 0.035 d -1 (at 235 0°C) for freshly produced detritus in the upper mixed layer. This "k" decreases with depth as an emergent result of the RC parameterization. The flux feeding rate depends on the particles' sinking flux, and thus attenuates the flux of large particles more strongly. Additionally, a small fraction of the large detritus intercepted by flux feeders is fragmented into small detritus. Phytoplankton growth rates and remineralization rates increase with temperature with a Q 10 of 1.9, whereas zooplankton growth rates have a Q 10 of 2.14. 240 To evaluate PISCES simulations against in situ POC measurements or their proxies, the correspondence between PISCES tracers (in italics) and observed POC fractions has to be established. In this study we assumed that sPOC corresponds to the sum of PISCES-simulated nanophytoplankton (PHY), diatoms (PHY2), small detrital particles (POC) and microzooplankton (ZOO), whereas bPOC corresponds to the sum of PISCES-simulated big detrital particles (GOC) and mesozooplankton 245 (ZOO2) ( Table 1). Total POC (tPOC) corresponds to the sum of those six PISCES tracers or, what is the same, sPOC + bPOC. Heterotrophic prokaryotes (BACT) are not a prognostic tracer in PISCES, and are not explicitly included in our analysis. The correspondence between observed and simulated POC fractions, explicit and implicit, is discussed in section 4.3.

PISCES 3D simulations vs. biome-aggregated observations 250
For the global-scale comparison between PISCES outputs and observations from BGC-Argo and satellites, we used the simulation presented by Aumont et al. (2017), which was forced by pre-computed dynamical fields from a pre-industrial run comparison with modeled fields, BGC-Argo observations were further screened to remove "outliers" in each biome and season. Outliers were defined as grid cells where the mean b bp700 in the upper 50 m was above the 95% percentile or greater than 0.008 m -1 (Briggs et al., 2020). The same resampling was applied to satellite-retrieved POC.

PISCES 1D simulations vs. BGC-Argo coherent annual time series (CATS)
A more detailed comparison was undertaken by matching each of the CATS from individual BGC-Argo floats with a 260 PISCES water-column ("PISCES 1D") simulation. The match was based on the coherence between the seasonal cycle of MLD observed by the float and the turbulent layer simulated by NEMO. The pre-computed dynamical fields used to evaluate the match-ups were obtained from an ocean-only historical simulation (NEMO v3.6) at 1-degree resolution with 75 vertical levels (ORCA1_L75 grid) forced with the JRA-55 atmospheric reanalysis that covered the 1958-2018 period (Tsujino et al., 2020) following the OMIP2 protocol (Griffies et al., 2016). The float-observed MLD and the NEMO-265 simulated turbocline depth (defined by turbulent diffusivity > 5·10 -4 m 2 s -1 ) were compared over the entire annual cycle in all the model grid cells that had been visited by the float on a given year. The best model grid cell was selected on the basis of model-observations correlation, root-mean-square error, and time lag in the onset date of summer stratification. An example of the metrics used to match NEMO dynamical fields and BGC-Argo MLD is provided in Fig. S1-S4.

270
PISCES 1D offline simulations were forced using the dynamical fields from the selected model grid cells. The same annual forcing, corresponding to the year of the BGC-Argo observations, was repeated over 5 simulation years. After 4 years of spin-up, the output from year 5 at 5-days resolution was used for the comparison to the BGC-Argo CATS. Initial conditions (climatological fields of inorganic nutrients and carbon chemistry variables) and boundary conditions (atmospheric deposition) were the same used for PISCES 3D (Aumont et al., 2015(Aumont et al., , 2017. Nutrient fields were restored towards the mean 275 annual profile below 300 m. This procedure avoided drift in nutrient stocks by replenishing the upper ocean with the same amount of nutrients each year, resulting in regular seasonal cycles after one year and identical cycles from year 4 onwards.

Climatological POC fields
This section describes the comparison among tPOC fields estimated from BGC-Argo and satellite observations and PISCES 280 simulations across four ocean biomes.

Permanently and seasonally stratified subtropical biomes
In the oligotrophic biomes, monthly median surface tPOC displayed low seasonal amplitude. Total POC concentrations 300 estimated from BGC-Argo data typically oscillated around 2 mmol m -3 , with a maximum/minimum ratio of around 1.6 in the Mediterranean and 1.3 in the Atlantic and Pacific subtropical gyres (STG). In the STG, satellite and BGC-Argo tPOC were in good agreement, which could be expected because b bp700 -POC conversion factors and satellite POC are based on the same study . On the contrary, PISCES tPOC exceeded BGC-Argo tPOC by around twofold in the STG. In the Mediterranean, satellite tPOC exceeded BGC-Argo tPOC by ~80%, pointing to the differences in the respective POC 305 estimation algorithms. PISCES tPOC was nearly fourfold higher than BGC-Argo tPOC at the surface in the Mediterranean, an overestimation that results from unrealistic physics in that basin caused by the coarse (2°) model grid (see 3.3 and 4.3).
Vertical tPOC profiles evidenced the shortcomings of PISCES simulations in the oligotrophic gyres. Compared to BGC-Argo profiles, PISCES simulations produced too-sharp deep POC maxima and underestimated tPOC in the waters above and below (Fig. 4). Overall, these patterns prompted us to examine in greater detail the seasonal cycles of POC in different 310 biomes.

Coherent annual time series of sPOC and bPOC: case studies
In this section we describe two CATS from BGC-Argo floats and their PISCES 1D counterparts (section 3. epipelagic tPOC all year round, and highest bPOC fractions of nearly 20% were recorded at the apex of the spring-summer blooms. Distinct vertical particle export events were observed in May and June, matching the surface phytoplankton blooms, and August, when nutrient limitation likely triggered bloom collapse. These export pulses produced synchronous increases in sPOC and bPOC through the mesopelagic layer, though with different magnitudes. After reaching relative minima in 335 October, mesopelagic sPOC and tPOC increased again in November, but they showed different vertical patterns. The matching PISCES 1D simulation captured with good skill the patterns of sPOC and bPOC in the epipelagic and, to a lesser extent, the mesopelagic layer (Fig. 6). Excellent correlation (r = 0.92) and bias (-0.5%) between BGC-Argo and PISCES were found for vertically integrated epipelagic tPOC (Fig. 6i). A delayed start of the bloom was observed in the 340 simulation, which can be partly attributed to a delay of around one week in the onset of permanent stratification in the model. It is also plausible that modelled phytoplankton reacted too weakly to the cessation of deep convection, which was captured by alternative MLD metrics in BGC-Argo profiles ( Fig. 6a and b). Despite the general good agreement between observed and simulated sPOC, the model produced a conspicuous plume of sPOC that sank from the surface spring bloom into the mesopelagic layer, at the prescribed constant rate of 2 m d -1 , which was not found in the observations. In the core of 345 this plume, PISCES POC exceeded BGC-Argo sPOC by more than twofold. A similar model-observations mismatch was observed in all the northern and southern subpolar CATS as well as in some CATS in the Mediterranean (Figures S5-S8 and https://doi.org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c Author(s) 2021. CC BY 4.0 License. S10). On average, simulated bPOC exceeded BGC-Argo bPOC by 36% and 96% in the epi-and mesopelagic layers, respectively. The largest overestimation was observed during the midsummer export event. On the other hand, bPOC underestimation was found during the May bloom between 0-400 m. 350 PISCES qualitatively reproduced the late summer peak of mesopelagic bPOC, which was observed in all the subpolar North Atlantic CATS. On the contrary, it failed to reproduce both the decrease of sPOC and bPOC in fall between 600-800m and the bPOC increase below 800 m. The latter occurred in 6 out of 11 CATS in the NASPG, all located in the Labrador Sea.
The apparent decoupling of deep mesopelagic bPOC from the overlying water column may be related to the insufficient 355 temporal resolution of BGC-Argo profiling during that period compared to bPOC sinking speed, or reflect bPOC export events from surface waters areas not located vertically over the float (Siegel and Deuser, 1997).

South Pacific subtropical gyre
Float 6901660 was deployed in March 2015 in the western South Pacific STG, and drifted westwards until it deflected SW while approaching Tahiti. As of March 2021, the float was still active and had completed 244 cycles with stable continuous 360 records at the 1000 m drift depth. Between July 2017 and June 2018 (18-21°S latitude and 148-157°E longitude), the period selected for the CATS analysis, the BGC-Argo profiles portrayed a stably stratified water column typical of the core of the subtropical gyres (Fig. 7). Vertical mixing events that reached a depth of around 100 m were observed in July, August and October. However, their effect on surface sPOC was hardly noticeable, indicating that turbulent entrainment of nutrients was too weak to stimulate new production significantly. A deep Chla maximum was present all year round between 150 and 200 365 m as identified by the maximum Chla gradient. This Chla maximum did not translate into a deep POC maximum. The fraction bPOC/tPOC was consistently around 6% in the epipelagic and 12% in the mesopelagic according to the vertically integrated stocks. The vertical-temporal distribution of bPOC was patchy in the lower mesopelagic (500-1000 m), perhaps reflecting the difficulty to detect rare aggregates from b bp700 spikes.

370
The epipelagic tPOC stock (driven by sPOC) simulated by PISCES matched well the observations, with a high temporal correlation coefficient (r = 0.69; Fig. 7i) despite the low seasonal variability in this tropical setting. The low mean bias of epipelagic sPOC (9%) resulted from mutually compensating biases, as PISCES overestimated sPOC between the base of the mixed layer and the deep Chla maximum, and underestimated sPOC below it. In the mesopelagic layer, vertically integrated PISCES sPOC was on average 45% lower than BGC-Argo sPOC, and showed low negative temporal correlation to the 375 observations. Regarding bPOC, the simulated stock was nearly twice as large as BGC-Argo estimates in the epipelagic, with the largest overestimation seen in the deep Chla maximum. In the mesopelagic, PISCES bPOC exceeded BGC-Argo estimates by 25% on average. Despite model-observations discrepancies in terms of bPOC concentration, PISCES simulations supported the increase in the bPOC/tPOC fraction between the epipelagic and the mesopelagic layers.

Coherent annual time series of sPOC and bPOC: generalized approach 380
Although each of the 50 CATS included in this study has unique features, some of the misfit patterns between PISCES and BGC-Argo data described in the previous section are common to the majority of CATS in a given biome or, more broadly, in subpolar vs. subtropical biomes. In this section we generalize the quantitative comparison between the 50 BGC-Argo CATS and their PISCES 1D counterparts across the four biomes. We consider separately the epipelagic and mesopelagic domains, focusing on mean annual sPOC and bPOC standing stocks (Fig. 8), the sPOC/tPOC fraction (Fig. 9), and the mesopelagic 385 Teff of sPOC and bPOC (Fig. 10).
Epipelagic sPOC stocks ranged between 193-425 mmol C m -2 and 282-537 mmol C m -2 , respectively, according to BGC-Argo observations and PISCES 1D simulations. Although the ranges of the different biomes overlapped, smaller stocks were usually found in the more oligotrophic (STG and Mediterranean) biomes. The agreement between simulated and observed 390 epipelagic sPOC was better in the STG and NASPG biomes, whereas simulated epipelagic sPOC exceeded observations by around 50% in the Mediterranean and Subantarctic biomes. Still, a significantly positive correlation between simulated and observed stocks (r = 0.45, p = 9·10 -4 ) was found overall.
In the mesopelagic domain, observed and simulated sPOC stocks ranged between 145-283 mmol C m -2 and 105-308 mmol 395 C m -2 , respectively, and the correlation between simulated and observed mesopelagic sPOC stocks was not significant (r = -0.08, p = 0.55). PISCES sPOC was similar to or slightly higher than BGC-Argo sPOC in the subpolar biomes, but up to twofold lower in the STG and Mediterranean biomes (as already shown in Fig. 4 and 7).

435
Our observational estimates are sensitive to the choice of b bp700 -POC conversion factors. These factors are especially uncertain in the mesopelagic due to data scarcity, which prompted us to use a constant global value for "c" in Eq 1. (Fig. 2); c = 1000 mmol C m -2 is the asymptotic POC-b bp700 conversion factor below 1000 m, taken from Cetinić et al. (2012). To address this uncertainty, we conducted sensitivity tests where c was halved or doubled. The resulting range of 500-2000 mmol C m -2 is probably generous, as our indirect estimates based on the study of Bishop et al. (1999) suggest a plausible 440 range of 815-1630 mmol C m -2 (Appendix A). As expected, changing c had little effect on epipelagic POC, a larger effect on mesopelagic POC and Teff, and no effect on the sPOC/tPOC fraction (because the same conversion factor is used to estimate sPOC and bPOC). Halving c resulted in steeper POC profiles, which brought mesopelagic sPOC closer to the 1:1 line in the STG and Mediterranean, at the expense of increasing the sPOC bias in the subpolar biomes and that of bPOC everywhere (Fig. S11). Doubling c caused a less steep vertical decrease of BGC-Argo POC, which overall worsened the model-445 observations agreement in the mesopelagic (Fig. S12) for sPOC stocks, bPOC stocks (except in the NASPG) and Teff.

Towards a globally consistent picture of POC fields in observations and models
Global quantification of POC stocks through the water column has been elusive until recently. On the one hand, direct chemical measurements (Lam et al., 2011) and indirect bio-optical observations from ships (Bishop et al., 1999;Gardner et 450 al., 2006) were sparse, and the ocean interior was severely undersampled compared to the euphotic layer. On the other hand, biogeochemical models were unable to capture even the order-of-magnitude POC concentration in the meso-and bathypelagic layers (Aumont et al., 2017). Our joint analysis of PISCES simulations, satellite observations and over 70,000 BGC-Argo vertical profiles reveals a globally consistent picture across the epi-and mesopelagic layers (Fig. 3-10 Table 2). The PISCES-derived estimate for the epipelagic layer, 1.6 Pg C, lies well within the range of previous estimates based on satellite observations. Gardner et al. (2006) estimated the global POC stock within the first attenuation depth using a compilation of in-situ POC and c p measurements leveraged by ocean color satellite data. They obtained a global stock of 0.43 Pg C, and invoked some scaling arguments to estimate that the total 465 POC stock down to middle-mesopelagic "background levels" would range 1-2 Pg C. Stramska (2009) obtained a global epipelagic POC stock of 1.8-2.3 Pg C using the algorithm of Stramski et al. (2008). Our analysis of the match-ups between BGC-Argo and satellite observations indicates that the latter algorithm overestimates POC at high latitudes outside of the summer season (Fig. 3). Potential explanations are the satellite algorithm being calibrated mostly against samples from lower latitudes (50°N-30°S), or its sensitivity to the differential atmospheric attenuation of the blue and green wavelengths at low 470 solar elevation (i.e. at high latitudes in winter). Recently, Evers-King et al. (2017) calculated the global mixed-layer POC stock using several satellite algorithms (including that of Stramski et al., 2008) and indicated that 0.77-1.3 Pg C was a plausible range. Our PISCES-based estimate for the mixed-layer POC stock, using the same MLD climatology (Schmidtko et al., 2013), is 0.58 Pg C. In this case, the lower PISCES-derived estimate may arise from the combination of POC overestimation by some satellite algorithms, discussed above, with PISCES' tendency to underestimate mixed-layer POC in 475 oligotrophic areas. This negative bias in PISCES is compensated by POC overestimation in deep Chla maxima below the MLD ( Fig. 4 and 7), resulting in good agreement between BGC Argo and PISCES over the entire epipelagic layer in the STG biome ( Fig. 8 and S13). Given the limited coverage of in situ seawater sampling and satellite observations in some regions and seasons (Evers-King et al., 2017), further intercomparison between those observations, BGC-Argo data and models is needed to better constrain epipelagic POC stocks. 480 POC in the lower mesopelagic and below has been traditionally treated as a "background" signal (Bellacicco et al., 2019;Gardner et al., 2006;Loisel and Morel, 1998). This approach is convenient for studies that focus solely on upper-ocean processes because POC concentration decreases exponentially with depth. Yet, our global estimates and several previous studies highlight the need to turn our attention to the large POC stocks (>2 Pg C) that reside in the mesopelagic and below it, 485 whose dynamics are still poorly understood. In line with previous studies (Dall'Olmo and Mork, 2014;Poteau et al., 2017), we showed that mesopelagic POC exhibits clear seasonal cycles in productive regions (Fig. 6, S5-S7 and S10), owing to their connection with the upper-ocean through numerous biological and physical processes (Boyd et al., 2019;Briggs et al., 2020). Despite being less reactive on average than upper-ocean POC, meso-and bathypelagic organic particles are microbial hotspots that host key biogeochemical functions, from enzymatic decomposition of macromolecules (Arnosti et al., 2012;490 Baltar et al., 2010a490 Baltar et al., , 2010b to aerobic and anaerobic respiration (e.g., Karthäuser et al., 2021) and chemosynthesis Pachiadaki et al., 2017). Moreover, mesopelagic particles are consumed by upper trophic levels that sustain fisheries (Bode et al., 2021;Woodstock et al., 2021).
In situ bio-optical measurements are poised to play a key role in monitoring marine POC stocks in layers that cannot be 495 accessed by remote sensing. For example, Sauzède et al. (2020) illustrated the power of merging BGC-Argo and satellite observations to obtain a dynamic 3D view of particle backscattering. Using a data-driven machine learning approach, they were able to predict the profiles of log 10 b bp700 measured by two BGC-Argo floats in the NASPG and STG biomes (R 2 of around 0.85 and MAPD of around 12%) from the sole knowledge of physical properties of the water column and surface ocean color (remote sensing reflectance). Their approach was recently extended to provide POC estimates (Sauzède et al., 500 2021), which can be of great utility for constraining biogeochemical models. Here we took an entirely different approach, based on converting b bp700 to POC with a simple empirical algorithm (Fig. 2) and then comparing it to the outputs of the PISCES model. Our PISCES-based estimates obtained a median R 2 = 0.86 and MAPD = 38% for 5-day depth-binned log 10 POC for 28 globally distributed floats. This good skill is remarkable because neither the empirical POC estimates nor PISCES were tuned to maximize their mutual agreement. Still, our study shows that the comparison of bio-optics-derived 505 POC measurements and PISCES is affected by different types of uncertainty that we analyze in the four sections below.

Bio-optical underpinnings of POC fields based on BGC-Argo observations
Accurate interconversion between bio-optical variables and concentrations is key for constraining ocean particle dynamics and their model representation (Bishop et al., 2004;Gardner et al., 2006). The variability in b bp700 -POC conversion factors is shaped by several concurrent processes that alter particle abundance, size distribution, shape, composition, and ultimately 510 https://doi.org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c Author(s) 2021. CC BY 4.0 License. optical properties. Below we discuss the processes that appear to drive, to first order, the variability of the POC/b bp700 ratios embodied in Eq. 1 and 2 ( Fig. 2; Appendix A), and the main strengths and weaknesses of our scheme.
Changes in the trophic status appear as the primary driver POC/b bp700 variability in the epipelagic layer (Cetinić et al., 2012; Figure A1). Productive waters host greater absolute and relative abundance of diatoms (Uitz et al., 2006) (see also  (Menden-Deuer and Lessard, 2000) and are covered with silica frustules that may scatter light more efficiently than naked cells (Twardowski et al., 2001), altogether resulting in lower POC content per unit b bp700 (Cetinić et al., 2012;Oubelkheir et al., 2005). The proportions of different autotrophic and heterotrophic organisms and detritus are also likely to vary with upper ocean productivity (see 4.3). If the mass-specific backscattering coefficients of these components were better known, their systematic variation patterns could be used to develop a continuous formulation 520 for POC/b bp700 , rather than the regionalized conversion factors used here. However, POC/b bp700 is influenced by other seawater constituents whose occurrence is less predictable. Foremost, biogenic calcite (e.g. from coccolithophores) and desert dust (Claustre, 2002;Loisel et al., 2011), both of which enhance b bp700 . In our data set, we detected a CATS (float WMO 6901647, year 2016) that was strongly affected by coccolith backscattering in the Iceland Basin, an area known for its massive coccolithophore blooms . This CATS was an obvious outlier in our model-observations 525 scatterplots (Fig. S13), likely because coccolith-enhanced b bp700 resulted in spuriously high observed tPOC, and was therefore removed from the analysis.
The decrease in POC/b bp700 along the vertical axis (Fig. 2) reflects the increase in the particles' index of refraction, hence the backscattering ratio (Cetinić et al., 2012;Nencioli et al., 2010). This change is likely caused by the remineralization of 530 organic materials (Martin et al., 1987) that leaves a higher mineral fraction (Honjo et al., 2008;Lam et al., 2011), and the increase in the structural complexity of aggregates with depth (Organelli et al., 2018). According to our sensitivity analysis ( Fig. S11 and S12), the prescribed exponential decrease of POC/b bp700 towards a constant POC/b bp700 (c = 1000 mmol C m -2 ) at depths >1000 m provides a good compromise globally, given the limited knowledge of mesopelagic POC/b bp700 and its variability across regions. However, the comparison between simulated and observed mesopelagic sPOC 535 (and thus tPOC) is more favorable in subpolar than in subtropical biomes. In the latter, better model-observations agreement was found when POC/b bp700 was halved (c = 500 mmol C m -2 ). A lower c would also bring the simulated Teff for sPOC and bPOC closer to BGC-Argo estimates. It is tempting to hypothesize that lower latitudes have lower mesopelagic POC/b bp700 owing to the greater proportion of calcite (Francois et al., 2002;Honjo et al., 2008;Lam et al., 2011Lam et al., , 2015. These hypotheses need further verification. 540 Finally, the decrease of surface POC/b bp700 with deeper vertical mixing imposed by Eq. 2 reflects the dilution of surface particle assemblages by entrainment of deeper waters (Lacour et al., 2019) with lower POC/b bp700 . wide seasonal MLD amplitude such as the NASPG. Note also that if POC/b bp700 in the mixed layer was constant, it would 545 have to decrease abruptly below the mixed layer to meet the low POC/b bp700 in the mesopelagic. On the other hand, the behavior prescribed by Eq. 2 may not be appropriate for situations when vertical mixing cannot erode the seasonal thermocline or pycnocline, because in such cases it will entrain water from the deep Chla maximum, and not mesopelagic water, to the upper mixed layer.

550
Our POC-b bp700 conversion algorithm omits several aspects that may have caused further uncertainty. Foremost, it assumes that sPOC and bPOC have the same POC/b bp700 ratio, which largely corresponds to that of the more abundant sPOC (Fig. 6-8). In addition, the in-situ data used to parameterize Eq. 1 and 2 are not free of uncertainty (Cetinić et al., 2012;Lam et al., 2011;Morán et al., 1999;Organelli et al., 2018;Strubinger Sandoval et al., 2021). While the development of more sophisticated POC-b bp700 interconversion schemes is desirable, the POC-b bp700 conversion scheme used here provides POC 555 estimates that are generally consistent with PISCES ( Fig. 8) and with previous assessments (Aumont et al., 2017). Moreover, this scheme allowed us to identify potential shortcomings of satellite-based assessments (Evers-King et al., 2017;Stramski et al., 2008). Indeed, bio-optical estimates of POC, and more generally BCP studies, would greatly benefit from the measurement of mass-specific bio-optical properties of various seawater constituents across different ocean basins and depths. 560

Correspondence between observed and simulated POC fractions
An additional source of uncertainty in our analysis is the imperfect match between the PISCES tracers and the observable POC fractions (Table 1). A positive aspect of our evaluation is the reasonable agreement between the simulated sPOC/tPOC fraction and the BGC-Argo estimates based on the b bp700 spike signal (Fig. 9). Our estimates converge to a global median value of 85%, with 95% of the data between 69-92%, fully within the range of size-fractionated chemical POC 565 determination (Aumont et al., 2017 and references therein). Our sPOC/tPOC estimates can also be compared to the suspended POC estimated with the marine snow catcher during spring in subpolar epipelagic and upper mesopelagic waters (Baker et al., 2017). The slow-sinking POC measured by the marine snow catcher should not be included in this comparison because it sinks at around 18 m d -1 , a range more typical of particles > 100 µm (Giering et al., 2016). The median 94% of suspended (non-sinking) POC reported by Baker et al. is higher than the 86% (81%) sPOC/tPOC estimated here from BGC-570 Argo (PISCES) in the subpolar biome, suggesting further comparison between different approaches is needed. The tendency of sPOC/tPOC to decrease between the epipelagic and the mesopelagic, found in both PISCES and BGC-Argo data (Fig. 9, Table 2), was also reported in previous studies (Aumont et al., 2017 and references therein;Organelli et al., 2020), lending further confidence to our estimates. On the other hand, the simulated partitioning of sPOC and bPOC into different living and detrital compartments is probably less realistic, as discussed in greater detail below.
A problematic aspect in the current representation of sPOC in PISCES is the partitioning between phytoplankton and nonphytoplanktonic (detrital + heterotrophic) POC in the upper ocean. This partitioning is far from being well established, neither from bio-optical observations nor from direct determination. Starting with bio-optical approaches, Oubelkheir et al. (2005) found that detritus accounted for around 60% of POC across a wide range in ocean productivity, and 70-85% of POC 580 was assigned to detritus + heterotrophs. These percentages accord well with those found by Claustre et al. (1999) in the tropical Pacific (their Table 1; see also Organelli et al., 2020). By contrast, Bellacicco et al. (2020) suggested that the nonphytoplanktonic b bp700 fraction varies widely, between <10% in the productive NASPG to >80% in subtropical gyres. Given the similar refractive indexes of living cells and organic detritus, hence similar POC/b bp700 , the latter and former estimates are hardly compatible when compared in POC units. Estimates based on direct POC stock determinations across wide 585 productivity gradients are also conflicting. Gasol et al. (1997) found a general increase in the autotrophic/heterotrophic ratio with productivity, with ratios as low as 0.1 at low autotrophic biomass. By contrast, Graff et al. (2015) found that phytoplankton accounted for between 15-85% of POC, with decreasing phytoplankton contributions at lower POC or Chla concentration (their figures 4 and 7), leaving no room for heterotrophs and detritus. Indeed, the two latter studies agree in the wide scatter in large-scale POC partitioning patterns, may further confound in situ bio-optical estimates. The low fraction of 590 detrital POC in PISCES near the sea surface, 15-31% ( Fig. 4; Table 2) is at odds with most observations, and suggests too slow transfer of planktonic biomass to detritus and/or too fast detrital POC removal, especially in oligotrophic waters.
PISCES bias may also arise from inappropriate representation of some POC reservoirs, such as heterotrophic prokaryotes (bacteria and archaea, BACT in PISCES), which have long been recognized as important contributors to the suspended POC 595 (Morel and Ahn, 1990;Gasol et al., 1997). However, BACT are not explicitly modelled in PISCES as prognostic tracers, meaning they are not interacting fully with other tracers. Instead, they are diagnosed in the productive surface layer from zooplankton biomass, based on an old model version that had interactive BACT. Below this layer, BACT biomass is propagated downwards with a power function based on Arístegui et al. (2009), which resembles a Martin curve (Martin et al., 1987) and is therefore very sensitive to the reference depth (Z 0 ) (Buesseler et al., 2020). We find two main arguments 600 against the inclusion of PISCES-estimated BACT in our POC estimates with the current model configuration. First, the empirical BACT estimation in PISCES has not been validated, to our knowledge, and may introduce noise in the comparisons. Second, PISCES-simulated POC already includes heterotrophic prokaryotes, because their biomass was not removed from the in situ POC measurements used to adjust the POC parameters in PISCESv2 (Aumont et al., 2017). In consequence, adding BACT to sPOC causes overestimation of mesopelagic POC (Fig. S14-S16) and can produce unrealistic 605 temporal patterns (Fig. S15). Nevertheless, we believe that inclusion of prognostic bacteria would enable more realistic simulation of POC stocks, with the potential side effect of improving the simulation of element fluxes in PISCES.
The comparison between simulated and observed bPOC is a novel contribution of our study. The method of Briggs et al. (2011Briggs et al. ( , 2020, originally developed to study intense POC export events, was here extended to estimate sPOC and bPOC separately over the full annual cycle through the epi-and mesopelagic domains. Our analysis shows that, despite the mismatch in terms of concentration, the bPOC derived from the spikes of high-resolution bio-optical profiles is strongly correlated to the PISCES-simulated bPOC (r = 0.78, r 2 = 0.61) along a wide trophic gradient (Fig. 8). This result is encouraging and supports the more widespread deployment of instruments that perform high-resolution bio-optical sampling to shed light on the spatiotemporal dynamics of large aggregates and particles (Briggs et al., 2020;Lampitt et al., 1993;615 Stemmann et al., 2008). On the other hand, it is unclear to what extent the bPOC inferred from the b bp700 spike signal is capturing mesozooplankton biomass, in addition to aggregates. Exclusion of PISCES mesozooplankton (ZOO2) from the comparison increases model-observations mismatch, with BGC-Argo bPOC exceeding PISCES estimates by around twofold, although the correlation remains nearly unchanged (r = 0.76, p < 10 -10 ). Imaging devices mounted on BGC-Argo floats may provide more accurate quantification of bPOC, allowing for the separation of detrital bPOC (Trudnowska et al., 620 2021) from mesozooplankton and micronekton (Haëntjens et al., 2020).

Importance of realistic physics and model evaluation across scales
In Fig. 5, PISCES simulations and BGC-Argo observations are compared using an array of skill metrics computed on the seasonal vertical profiles of tPOC between 0-1000 m. Starting with the 3D seasonal climatology, we observed that the correlation between the median (aggregated) profiles was generally better than the correlation between the spatially 625 collocated (non-aggregated) profiles within a given biome and season (Fig. 5a), whereas no differences were found in terms of dispersion metrics (Fig. 5b and c). The difference in correlation was larger in subpolar biomes, suggesting that the modelobservations spatial mismatch was magnified in regions with more energetic ocean dynamics and sharper physical and biogeochemical gradients, whose real-world location may not be well reproduced by the ocean circulation model used to force PISCES. For PISCES 1D simulations, their correlation coefficients with their CATS counterparts was usually in the 630 high range of the correlation coefficients obtained by the biome-median PISCES 3D profiles. In terms of dispersion metrics, the ensemble of PISCES 1D simulations showed wider dispersion, but the best 1D simulations clearly outperformed the 3D simulation in a given region and season. The better skill of 1D simulations was more evident during spring, a season characterized by the onset of stable stratification after deep winter vertical mixing in middle and high latitudes. The greatest difference between 3D and 1D simulations was found in the Mediterranean, highlighting the more realistic vertical mixing 635 and upper ocean productivity in the 1D simulations.
Our cross-scale evaluation indicates how crucial it is to evaluate model physics before extracting conclusions on biogeochemical model performance (Doney et al., 2004;Kriest et al., 2020;Löptien and Dietze, 2019). In our 1D CATS approach, the skill of PISCES simulations was maximized by carefully matching observed and modelled vertical mixing 640 ( Fig. S1-S4), which is a key driver of upper ocean ecosystems. This approach has a subjective component, and may also suffer from the idealized assumption that BGC-Argo profiles reflect mostly vertical-scale processes, disregarding horizontal advection (Alonso-González et al., 2009). Yet, the similar misfit patterns encountered for different CATS within a given biome support the robustness of the 1D matching approach (compare Fig. 6 versus S6, and Fig. 7 versus S9). To further evaluate this issue, we matched different neighboring NEMO grid cells to the same in situ CATS. Again, this exercise 645 indicated that our main conclusions are not sensitive to the choices made for 1D model-observations matching (compare Fig.   6 and Fig. S5). Indeed, alternative matching approaches can be devised, each of which with advantages and pitfalls, for example: (1) sampling the outputs of biogeochemical models at the locations visited by BGC-Argo floats, which may require high resolution models; (2) deploying virtual BGC-Argo floats (van Sebille et al., 2018) and comparing them statistically to observations; or (3) forcing 1D biogeochemical simulations with observed physical fields, e.g. vertical mixing (Llort et al., 650 2015) or light (Terzić et al., 2019). As a general rule, the good skill of the best PISCES 1D simulations (Fig. 5) indicates that our CATS approach can be used to tease apart model-observations misfits caused by model physics from those caused by the parametric uncertainty and the structure of the biogeochemical model, opening up new avenues for parameter optimization (Falls et al., submitted) and model development.

655
Our cross-scale evaluation is also informative as to the spatiotemporal scales that can be addressed with a given model setup.
This matter was recently tackled by Bisson et al. (2019) using the export production model of Siegel et al. (2014), which is forced by satellite-derived primary production. In particular, they showed that this diagnostic model could be optimized to reproduce global climatological patterns, but exhibited poor skill when faced with non-climatological datasets, which reflect local snapshots of ecosystem functioning. Model evaluation at climatological scales provides an incomplete picture, 660 especially in productive regions where much POC export can take place during intense but short-lived events (Briggs et al., 2020). Such events are smoothed out when climatologies are computed, and their coherence with the physical forcing is lost.
Our work shows that a prognostic model like PISCES can afford both event-scale and climatological scale predictions. This capability is important to test our process-level understanding, which underpins climate change projections.

Joint use of BGC-Argo and models for process-level understanding 665
Properly representing POC stocks is crucial for constraining epi-and mesopelagic carbon budgets and, ultimately, estimating the strength of the BCP and predicting its future evolution. The mismatch patterns between simulated and observed POC profiles (Fig. 4, 6 and 7) indicate different types of model shortcomings in subpolar and subtropical latitudes. The poorest agreement between PISCES and BGC-Argo data is found when their respective estimates of mesopelagic Teff are compared ( Fig. 10), which should prompt further research on this key descriptor of the BCP, both in terms of POC fluxes (Buesseler et 670 al., 2020) and stocks (Lam et al., 2011;this study), and their relationship with the structure and productivity of the upper ocean ecosystem.
In the subpolar biomes, we find discrepancies between the patterns of sPOC and bPOC vertical export triggered by the intense surface phytoplankton blooms typical of these waters (Fig. 6, S5-S7 and S10). The observed rapid sPOC increase at 675 depth cannot be explained by the sPOC gravitational sinking, and fragmentation of rapidly sinking bPOC has to be invoked https://doi.org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c Author(s) 2021. CC BY 4.0 License. (Briggs et al., 2020). This fragmentation is most likely caused by zooplankton feeding Stemmann et al., 2004aStemmann et al., , 2004bStukel et al., 2019) and swimming (Goldthwait et al., 2004), combined with bacterial hydrolysis of aggregatebinding polymers (Arnosti et al., 2012;Baltar et al., 2010a). In relative terms, mesopelagic sPOC increases less strongly than bPOC during blooms (Fig. 6i), which is consistent with sPOC being a byproduct of the transformation of less-abundant 680 bPOC. Fragmentation processes may supply fresher sPOC to the mesopelagic, enhancing the coexistence of suspended particles with variable freshness (Alonso-Gonzalez et al., 2010;Aumont et al., 2017) and overall contributing to POC remineralization (Giering et al., 2014;Mayor et al., 2020).
In the subtropical gyres and the most oligotrophic Mediterranean waters, PISCES underestimates tPOC between 0-1000 m 685 except for the deep Chla maximum, where it overestimates tPOC. The prominent deep POC maximum simulated by PISCES is generally not found in observations from the STG biome ( Fig. 7 and S9), where deep Chla maxima generally reflect phytoplankton photoacclimation, not enhanced phytoplankton biomass (Cornec et al., 2021). Thus, PISCES possibly overestimates the productivity of deep Chla maxima globally, indirectly causing stronger nutrient limitation at the surface.
Between the deep Chla maximum and 200 m, POC pools decrease more steeply in PISCES than in BGC-Argo observations. 690 Between 200 and 700 m, by contrast, simulated sPOC and bPOC Teff are only 10% lower than observations, well within observation uncertainty. Thus, the mesopelagic sPOC deficit simulated by PISCES in the STG originates mostly through the insufficient vertical POC export at 200 m depth. The Mediterranean represents a different case, whereby the large disagreement between PISCES and BGC-Argo Teff may result from either poorly constrained POC/b bp700 , simulated POC dynamics, or both. Thus, this region may provide a good testbed for studying the role of the mineral fraction in the BCP, 695 including the controversial ballast hypothesis (Francois et al., 2002;Klaas and Archer, 2002;Passow, 2004).
In summary, the mismatch with observations suggests the need to improve the representation of sPOC-bPOC interconversion in PISCES. The size distribution of POC along the vertical axis is a key variable for constraining POC budgets, because it reflects the interplay between gravitational sinking, remineralization, trophic transfer and 3D dynamics 700 including horizontal POC advection (Alonso-González et al., 2009;Boyd et al., 2019). In many instances (Fig. 9), our analysis suggests that better model performance in the mesopelagic may be achieved by increasing the net transfer of bPOC to sPOC and the Teff of both fractions, which may require further balancing the interplay between various mechanisms of POC export and flux attenuation. The vertical model-observations mismatch patterns observed here emphasize that POC budgets have to be computed with the highest vertical resolution affordable, or otherwise an apparent POC budget balance 705 may result from compensating imbalances in different horizons (Giering et al., 2014;Marsay et al., 2015). The detailed level of information available from BGC-Argo floats may prove to be extremely valuable to help improve the POC schemes embedded in models such as PISCES.

Conclusions and outlook
In this study we compared globally distributed POC observations between 0-1000 m made by BGC-Argo floats to the 710 predictions made by the PISCES model (PISCESv2). A subset of BGC-Argo floats profiling at high vertical resolution enabled us to analyze small and big POC separately. The comparisons rely on a proposed new scheme for converting a biooptical measurement (b bp700 ) to POC. Although PISCES recreates with good skill the main features observed in subpolar and subtropical biomes, the comparison is still hampered by: (1) spatial and temporal variability in POC/b bp700 conversion factors, (2) mismatches in observed and simulated physics, and (3) imperfect correspondence between observed and 715 simulated POC fractions. Evaluation of these uncertainties allowed us to detect limitations of the biogeochemical model parameterizations, which may arise from both suboptimal model parameters and model structure. Therefore, the descriptive work and model-data matching strategies presented here pave the way towards the use of BGC-Argo observations for data assimilation and model parameter optimization (Falls et al., submitted) and, ultimately, model development.

720
Widespread use of BGC-Argo data for understanding POC budgets and the BCP can complement classical model constraints based on vertical POC fluxes and ocean-interior nutrient remineralization. Merging of BGC-Argo and satellite data streams through data-driven approaches -which allow for great flexibility-, and mechanistic models -which provide processlevel understanding-, can soon provide us with a high-resolution 4D view of the oceanic carbon cycle. Further work is granted to investigate POC dynamics through combination of PISCES2 and autonomous observations. 725 Below we list several research priorities, whose implementation would advance the study of the biological carbon pump through the synergies between BGC-Argo observations and modelling. o Further developments in BGC-Argo data processing are needed, with the final goal of supplying a wide public with user-friendly data products in near real time.
• Models: o Evaluating model skill at resolving POC stocks, in addition to fluxes, is key to ensure that models 740 reproduce observed fluxes for the right reasons. o Continuous development of schemes representing particle dynamics across the entire size spectrum is needed to constrain ecologically, climatically and economically relevant element fluxes. 745 o Extension of prognostic modelling of bio-optical properties (Dutkiewicz et al., 2019) into the meso-and bathypelagic layers would enable direct matching with measurements made from autonomous platforms, facilitating their assimilation by models.
• Joint planning of field observation and modelling projects, from their very conception and through their entire development, is key to fully exploit the capabilities of each approach. 750 Data availability. The simulated and observed datasets analyzed in this article are available at https://doi.org/10.5281/zenodo.5139602. The code and documentation of NEMO and PISCES are available at 755 https://www.nemo-ocean.eu/. The authors can provide the code used to process the datasets on reasonable request.

Supplement.
A supplement with 16 figures (Fig. S1-S16) and one table (Table S1)  Appendix A: Calculation of POC/b bp700 ratios and related optical considerations POC/b bp700 ratios at the sea surface were obtained from the slope of the linear regression between POC and b bp700 (Table A1).
Our approach essentially follows the literature compilation made by Cetinić et al. (2012). The linear regressions generally yielded small and non-significant y intercepts for sea-surface data. Therefore, we assumed the slopes were equal to the 775 POC/b bp700 ratio. When the b bp measurements were made at another wavelength, we converted them to b bp700 assuming an exponent η = 0.41 (Cetinić et al., 2012), according to the equation: For measurements taken at 555 nm, which was the wavelength used in two studies (Table A1), this resulted in 10% lower b bp700 compared to b bp555 , thus higher POC/b bp700 ratios. 780 The dataset of Cetinić et al. (2012) was the only one showing a significant negative intercept of the POC vs. b bp700 linear regression. Unlike the other datasets, where the POC vs. b bp700 relationship reflected mostly sea-surface variability in a given biome, this was the only dataset that included data collected between the sea surface and 600 m. Bol et al. (2018) reprocessed this dataset by computing the linear regression between POC and b bp700 in different depth bins, forcing the 785 intercept through zero. The slope of the POC vs. b bp700 regressions, and therefore the POC/b bp700 ratio, decreased by more than twofold between the surface and the 600 m bins. These results suggest that the POC vs. b bp700 relationship may be better modelled with nonlinear equations, as done for surface data in some previous studies (Balch et al., 2010;Johnson et al., 2017;Stramski et al., 1999). However, one must keep in mind that the ecosystem processes that define POC/b bp700 ratios in the surface layer may be different from those occurring in the mesopelagic (see Discussion,4.2). 790 The conversion between POC and particulate beam attenuation at 660 nm, c p , in the surface ocean has been analyzed on more occasions than its b bp700 counterpart. The backscattering ratio, b bp /b p , relates particulate backscattering to the total particulate scattering, and is directly related to the refractive index of the particle assemblage (Babin et al., 2003;Dall'Olmo et al., 2009;Loisel et al., 2011;Organelli et al., 2018;Stramski and Kiefer, 1991;Stramski et al., 1999;Twardowski et al., 795 2001;Ulloa et al., 1994). Light absorption by particles is negligible in the 650-700 nm spectral region, such that total beam attenuation c p is a good approximation of the scattering coefficient b p . Thus, once the known absorption and scattering coefficients of seawater are removed, b bp /c p is a good approximation of the backscattering ratio, and can be used to compare the relationships between POC vs. c p and POC vs. b bp700 .

800
The observed POC/b bp700 variability in the surface layer across biomes reflected in our POC estimation algorithm is analogous to that found for the relationship between POC/c p by Cetinić et al. (2012), as shown in Figure A1, although fewer POC-b bp700 datasets are available. The POC content per unit b bp700 decreases with the maximum b bp700 of each dataset, which may be related to the structure and species composition of upper-ocean ecosystems (see Discussion,4.3). Given that the https://doi.org/10.5194/bg-2021-201 Preprint. Discussion started: 18 August 2021 c Author(s) 2021. CC BY 4.0 License.
backscattering coefficient is sensitive mostly to particles smaller than 30 µm (Organelli et al., 2018), it is also plausible that 805 the portion of POC that contributes to backscattering varies across sites, further enhancing POC/b bp700 variability.
The number of studies that tackled POC vs. bio-optical conversion in the mesopelagic layer is much smaller than those that focused on the epipelagic layer. Besides Cetinić et al. (2012) in the subpolar North Atlantic, we are only aware of the study of Bishop et al. (1999) in the northeast Subarctic Pacific. The latter study found a POC vs. c p slope of 16.3±0.49 mmol C m -2 810 (= mmol C m -3 / m -1 ) over the 0-1000 m depth range (the positive regression y-intercept in that study corresponds to absorption by seawater). This slope lies in the low range of the compilation shown in Fig. A1a (196±5.9 mg C m -2 in mass units). A close-up into the lower measurement range (Fig. 5 of Bishop et al., 1999) shows that the slope was not clearly distinct for POC concentration < 0.5 mmol m -3 , typical of the mesopelagic layer. If we assume that the backscattering ratio was likely between 1 and 2% in the mesopelagic, then we obtain a POC/b bp700 ratio of 815-1630. This range is compatible 815 with the POC/b bp700 of ~1000 mmol C m -2 found at 600 m by , which we used as the asymptotic deep POC/b bp700 ratio in Eq. 1. Note that the linear relationship between c p and POC is compatible with the nonlinear relationship between POC and b bp700 along the vertical profile if the backscattering ratio also increases with depth, as found by Cetinić et al. (2012). The latter study found an increase in the backscattering ratio (b bp700 /c p653 ) from around 1.2% at the surface to around 1.5% in the upper mesopelagic. The b bp700 /c p653 ratio was more variable in deeper layers and values > 2% were not 820 rare (see also Nencioli et al., 2010;Organelli et al., 2018;Twardowski et al., 2001).
Beyond the natural variability, the interconversion between POC and bio-optical proxies is also confounded by methodological, most of which were not fully addressed in the studies compiled here. Measurement of POC in the lower mesopelagic (below 500 m) requires filtration of large sample volumes in the m 3 range (Bishop et al., 1999), otherwise POC 825 concentration rapidly approaches the detection limit with procedures and filtration volumes used for the epipelagic and the upper mesopelagic (Cetinić et al., 2012;Strubinger-Sandoval et al., 2021).
In practice, it includes heterotrophic prokaryotes' biomass with current parameter values (see 4.3).
In practice, measurements may include particleattached and free-living organisms, viruses, colloids and adsorbed DOC.

Microzooplankton. Mostly ciliates and flagellates
with size similar to their prey, down to around 2 µm (Calbet 2008) Heterotrophic prokaryotes (BACT): currently not a prognostic tracer in PISCES. Not considered explicitly in this study (see 4.3).