BPOP-v1 model: exploring the impact of changes in the biological pump on the shelf sea and ocean nutrient and redox state

The ocean’s biological pump has changed over Earth history from one dominated by prokaryotes, to one involving a mixture of prokaryotes and eukaryotes with trophic structure. Changes in the biological pump are in turn hypothesised to have caused important changes in the ocean’s nutrient and redox properties. To explore these hypotheses, we present here a new box model including oxygen (O), phosphorus (P) and a dynamical 10 biological pump. Our Biological Pump, Oxygen and Phosphorus (BPOP) model accounts for two – small and large – organic matter species generated by production and coagulation, respectively. Export and burial of these particles are regulated by a remineralization length (zrem) scheme. We independently vary zrem of small and large particles in order to study how changes in sinking speeds and remineralization rates affect the major biogeochemical fluxes, and O and P ocean concentrations. Modelled O and P budgets and fluxes lie reasonably 15 close to present estimates for zrem in the range of currently measured values. Our results highlight that relatively small changes in zrem of the large particles can have important impacts on the O and P ocean availability and support the idea that an early ocean dominated by small particles was nutrient rich due to inefficient removal to sediments. The results also suggest that extremely low oxygen concentrations in the shelf can coexist with an oxygenated deep open ocean for realistic values of zrem, especially for large values of the small particle zrem. This 20 could challenge conventional interpretations that the Proterozoic deep ocean was anoxic, which are derived from shelf and slope sediment redox data. This simple and computationally inexpensive model is a promising tool to investigate the impact of changes in the organic matter sinking and remineralization rates as well as changes in physical processes coupled to the biological pump in a variety of case studies.


Introduction 25
The 'biological pump' describes the production of organic matter at the ocean's surface (an oxygen source), its downward export/sinking flux, remineralisation at depth (an oxygen sink), and burial. This set of processes acts against the homogenization of tracer concentrations by the ocean's circulation, maintaining large-scale tracer gradients (Sarmiento and Gruber, 2006). In today's world, the biological pump plays a key role in transferring carbon from the atmosphere/surface ocean to the deep ocean and in so doing lowers atmospheric CO2 and creates 30 oxygen demand in deeper waters (Lam et al., 2011;Kwon et al., 2009). Those deeper waters with the greatest oxygen demand relative to oxygen supply can be driven hypoxic (O2 < 60 mmol m -3 ), suboxic (O2 < 5 mmol m -3 ) or even anoxicas is being seen in parts of the ocean today (Keeling et al., 2010). By combining surface oxygen production and organic carbon burial, the biological pump plays a role in determining the long-term source of oxygen to the atmosphere. The biological pump also provides a means of efficiently transferring organic matter 35 and the nutrients it contains to marine sediments, if sinking through the water column happens fast enough compared to remineralization for the material to hit the bottom (Sarmiento and Gruber, 2006). Hence the biological pump plays a key part in balancing the input of phosphorus to the ocean with a corresponding output flux of phosphorus buried in marine sediments.
Through Earth's history, the characteristics, efficiency and impact of the biological pump are thought to have 5 changed dramatically due to the evolution of increasingly large and complex marine organisms (Ridgwell, 2011;Logan et al., 1995;Boyle et al., 2018). Life in the ocean began as just prokaryotes, presumably attacked by viruses, with slow sinking of the resulting tiny particles. Now the marine ecosystem is a mix of prokaryotic cyanobacteria and heterotrophs, and size-structured eukaryotic algae, mixotrophs and heterotrophs all the way up to large jellyfish, fish and whales. Some of the resulting particles sink very fast (McDonnell and Buesseler, 2010). 10 How changes in the biological pump have affected ocean nutrient and redox state at different times in Earth history is a subject of active research and hypothesis generation. Previous work has highlighted the Neoproterozoic Era, spanning from 1,000 to 541 million years ago, as of particular interest because it saw a shift of dominance from prokaryotes to eukaryotes and a series of dramatic shifts in the climate, biogeochemical cycling and ocean redox state (Katz et al., 2007;Brocks et al., 2017). A common paradigm has been to assume that a progressive rise of 15 oxygen in the atmosphere (of uncertain cause) drove the oxygenation of the deep ocean at this time through airsea gas exchange and mixing, but equally increases in the efficiency of the biological pump could have lowered ocean phosphorus concentration and thus oxygenated the ocean (Lenton et al., 2014). Recent data show a series of transient ocean oxygenation events ∼660-520 Ma, which get more frequent over time, suggesting a complex interplay of processes on multiple timescales, including changes in the biological pump and ocean phosphorus 20 inventory .
During the Phanerozoic Eon there have been further changes to the biological pump. In particular, a rise of eukaryotic algae from the early Jurassic onwards is hypothesised to have increased the efficiency of the biological pump and thus oxygenated shallow waters (Lu et al., 2018), but presumably deoxygenated deeper waters, at least in the short term. In the oceanic anoxic events (OAEs) that occurred during the Mesozoic Era there were major 25 increases in prokaryotic nitrogen fixation yet evidence for a eukaryote-dominated biological pump (Higgins et al., 2012), raising interesting questions as to whether this reinforced anoxia at depth.
Previous modelling work has examined the impact of changes in the organic matter remineralisation length/depth (zrem) in the 3D GENIE intermediate complexity model (Meyer et al., 2016;Lu et al., 2018). Both studies clearly demonstrated the important control of the zrem on ocean oxygen concentrationsas it gets larger the oxygen 30 minimum zone shifts to greater depths. Furthermore, Lu et al. (2018) showed that an increase in zrem can explain an observed deepening of the oxycline from the Paleozoic to Meso-Cenozoic in the ocean redox proxy I/Ca. However, coarse 3D models such as GENIE do not really resolve shelf seas and their dynamics, which are distinct from those of the open ocean. Furthermore, GENIE only accounts for one organic carbon species, overlooking processes of transformation of organic material, such as coagulation and fragmentation, which contribute to 35 modulate the efficiency of the organic matter vertical export and burial (Wilson et al., 2008;Karakaş et al., 2009;Boyd and Trull, 2007).
In this study we take a more idealised approach, exploring how changes in the properties of the biological pump may have affected the shelf sea and open ocean nutrient and redox state using a new Biological Pump, Oxygen and Phosphorus (BPOP) box model. This model combines a box representation of the marine O and P cycles with 40 an intermediate complexity representation of the biological pump transformations, including two classes of particulate organic matter (POM). BPOP allows us to modify the properties of two POM pools, whose abundance is regulated by the processes of production and coagulation. We focus on changes in the characteristic depths at which the two POM pools are remineralized, i.e., the particle remineralization length scale (zrem), and study the resulting equilibrium budgets and fluxes. In the following sections we describe the model, we provide an 5 evaluation of its performance in the context of modern observations and flux estimates, and finally present and discuss our model results.

Model description
Here we describe the Biological Pump, Oxygen and Phosphorus (BPOP) model. The model was implemented using Matlab and the equations are solved by the built-in ode15s solver. BPOP can easily run on a single core, 10 integrating 50 million years of time in less than a minute on an ordinary machine, and is therefore computationally efficient. We refer to the user's manual (see the supplementary material) for further information on how to run the model.

Variables and circulation
The box model resolves explicitly for each relevant box the local concentrations of three types of tracers: 15 molecular oxygen O2 (O), inorganic dissolved phosphorus (P) and sediment organic phosphorus (SedPorg). The total budgets of P and O, respectively P TOT and O TOT , are also independently integrated from the net sources and sinks of the two tracers over the entire model domain, for the purpose of checking mass conservation. The entire set of the model's state and diagnostic variables and their units are listed in Table 1. In the following subsections we describe the box model's geometry and discuss the physical and geochemical fluxes that drive the tracers ' 20 dynamics. Box properties are listed in Table 2, while the set of parameters adopted for the modelled physical and geochemical fluxes can be found in Table 3.

Box properties and physical fluxes of inorganic tracers
The box model includes 4 ocean boxes, 1 atmospheric box and 2 sediment boxes ( Figure 1a). The ocean and sediment boxes are equally split between shelf sea and open ocean, both including one surface ocean box and one 25 deep ocean box.
O and P are exchanged between the 4 ocean boxes through advection and mixing, including an upwelling recirculation between shelf sea and open ocean (Wollast, 1998). For a generic tracer concentration C and in the i th box, the physical exchange flux is represented by For each surface box i, air-sea gas exchange allows O fluxes between the ocean and the atmosphere (at). The flux is positive when directed into the ocean and depends on the gas transfer velocity Kw and Henry's constant KHenry, as in: where KW is a function of the prescribed mean temperature Tmean and wind speed Wspeed (Sarmiento and Gruber, 35 2006).

Initialization and boundary fluxes
The model is initialized with an even concentration of P (Pini) in all the ocean boxes, zero oxygen and zero sediment Porg. A constant input of P from rivers (Pin) into the surface ocean replenishes the P ocean reservoir despite the burial flux (net sink of Porg) into the sediments. Pin is in part delivered directly to the surface open ocean (Sharples et al., 2017). At equilibrium, the Porg burial flux balances Pin. Oxidative weathering determined 5 by atmospheric oxygen O at constitutes a net sink flux for O. The weathering flux depends on a constant a baseline flux W0 and it scales like the square root of the oxygen concentration normalised to present values Omix0 , following:

Biological pump details 10
The modelled tracer cycles are coupled by a set of biological transformations, i.e., the biological pump, governing the cycle of production, remineralization and burial of Porg in the water column and in the sediments. Porg in the water column is resolved implicitly: at each time step all the produced Porg that does not reach the sediments is instantaneously remineralized. In this sense, in our model no Porg can accumulate in the ocean's water column.
This scheme is like the one used to represent detrital POM in some modern ocean biogeochemical models (Moore 15 et al., 2004). P and O biological fluxes are coupled by a fixed Redfield ratio OPRed. The next few paragraphs describe the cycle of production, coagulation, export, remineralization and burial that constitute the biological pump representation. The full set of parameters used to resolve the Porg cycle is provided in Table 4.

Particle classes, production and coagulation
The model includes two Porg classes which get produced, exported and remineralized in the ocean's water column: 20 small Porg (SPorg) and large Porg (LPorg). The use of two Porg classes is in line with modern ocean in situ observations, which reveal a bimodal distribution of the particle sizes and sinking speeds (Riley et al., 2012;Alonso-González et al., 2010). Moreover, it allows to better reproduce the commonly observed Martin power-law decay of the particle export flux with the use of a remineralization length scheme of export and burial fluxes (Boyd and Trull, 2007). 25 Organic matter production happens only in the surface ocean boxes through the uptake of P. This is regulated by a maximum rate Peff and a Michaelis-Menten kinetics with constant KP. Production in each i th box only generates SPorg, according to: LPorg is generated via the coagulation of SPorg either at the surface after production or at depth after the export of 30 SPorg from the surface. The coagulation of SPorg into LPorg in each box i is proportional to the square of the local SPorg concentration and is regulated by a coagulation rate cgr (Boyd and Trull, 2007;Gruber et al., 2006), as in: Coagulation impacts the relative contribution of small and large particles to the export and burial fluxes by subtracting from the local SPorg pool and adding to the LPorg pool. 35

Physical fluxes of organic material
The implicit representation of the organic matter in the water column implies that no organic matter is accumulated in the ocean. In our baseline version of the model, corresponding to the results presented in this manuscript, SPorg and LPorg are redistributed throughout the watercolumn exclusively by implicitly modelled gravitational sinking before being either buried, accumulated in the sediments or remineralized. Even though the vertical export by 5 downwelling and mixing (Stukel and Ducklow, 2017), and the lateral organic matter redistribution (Lovecchio et al., 2017;Inthorn et al., 2006) may be important when working with suspended SPorg (zrem S = 0), these fluxes are not currently accounted for in the model.

Remineralization length scheme
The export and sedimentation fluxes of Porg through the water column are represented by a remineralization length 10 scheme. In this representation, the vertical fluxes of organic matter f(z) vary exponentially with depth. The shape of the exponential depends on the value of the remineralization length (zrem) of each organic matter species: where is the flux at the reference depth , and the index indicates the organic matter pool of reference, either small (S) or large (L). This representation of the export flux is convenient, as it does not depend on the specific 15 choice of (Boyd and Trull, 2007).
The remineralization length zrem indicates the distance through which the particle flux becomes 1/e times (about 36 %) the flux at the reference depth (Buesseler and Boyd, 2009;Marsay et al., 2015). This quantity is expressed in metres and can be calculated as the ratio between the particle sinking speed and the particle's remineralization rate . Consequently, zrem implicitly contains information on several particle inherent properties 20 (among which density, size, shape, organic matter liability) as well as information about the surrounding environment, e.g., the type of heterotrophs which feed upon the organic material (McDonnell and Buesseler, 2010;Baker et al., 2017). For simplicity, we assume that the remineralization length of small and large particles does not vary between shelf sea and open ocean boxes. We examine the potential impact of this limitation in the discussion section of the paper. 25

Sediments and burial
SPorg and LPorg accumulate in the sediments as SedPorg, which is calculated as a density per unit of area. The sediment flux into the sediment box i depends on the organic matter concentration in the overlaying deep ocean box j and on the remineralization length of the two pools as in: The accumulated SedPorg is partially slowly remineralized and partially irreversibly buried in a mineral form.
Phosphorus burial as mineral Ca-P is modelled as a function of the square of SedPorg that accumulates in the sediments, in a way that is analogous to the dynamics of particle coagulation, and is regulated by a constant rate coefficient CaPr. Ca-P formation happens at a lower rate under low oxygen conditions (CaPr * = CaPr • fsan with fsan < 1), in agreement with observations and previous models (Slomp and Van Cappellen, 2006). The transition 35 from aerobic and anaerobic conditions is controlled by a Michaelis-Menten type of function of the oxygen concentration in the deep ocean box j overlaying the sediment box i. The oxic and anoxic terms sum to the total formation term as in: This flux is essential to balance the continuous P river input, therefore preventing the ocean from overflowing with nutrients. 5

Remineralization in the water column and sediments
At each time step, remineralization in the water column completely depletes the Porg that has not reached the sediments. In the two surface boxes, remineralization of Porg that is not exported below the euphotic layer uses up part of the oxygen that was released by production. For this reason, net oxygen production in each surface box is proportional to the export of Porg below the euphotic layer. Export from a surface box i to a deep box j happens 10 both via gravitational sinking and via mixing, as in: At depth, the remineralization of Porg that does not reach the sediments happens through both aerobic and anaerobic processes, completely depleting the remaining Porg. Water-column remineralization of Porg into inorganic P in the deep box j is therefore calculated as: 15 In each deep ocean box i, aerobic remineralization uses some of the available oxygen and is therefore limited by a Michaelis-Menten kinetics with a half-saturation constant K w O (DeVries and Weber, 2017). Anaerobic remineralization takes up the entire remaining Porg that is not remineralized aerobically and releases a product which "bubbles up" to the atmosphere, reacting with atmospheric oxygen. 20 In each sediment box i, remineralization of SedPorg happens in a similar way to remineralization in the water column, with an aerobic and an anaerobic component. However, remineralization in the sediments is not instantaneous, but happens at a fixed rate which depends on the oxygenation state of the overlaying water column.
Aerobic remineralization takes up oxygen from the overlaying deep-water box j and happens at a rate rmr, while being limited by a Michaelis-Menten coefficient. Anaerobic remineralization releases its product to the 25 atmosphere and happens at a faster rate rmr * = rmr •fean with fean >1, in agreement with recent observations and previous models (Slomp and Van Cappellen, 2006). Total sediment remineralization is therefore the sum of the two terms as in:

Equations summary
The dynamics of the model's 11 state variables is regulated by just as many equations. We summarize here the major terms for P, O and SedPorg in the surface ocean (s), deep ocean (d), atmosphere (at) and sediments, without distinguishing between coastal and open ocean boxes and assuming that all terms have been scaled with the 35 reference box' dimensions. A full set of equations including the explicit formulation of all the flux terms for each box can be found in the paper's Appendix.

Strategy: sensitivity studies for varying zrem 15
In order to characterize the model, we analyse the equilibrium budgets and fluxes of the state variables for varying zrem values separately for SPorg and LPorg, respectively zrem S and zrem L . We adopt a range of zrem values that fall close to modern observations Buesseler and Boyd, 2009;Marsay et al., 2015) and keeps into consideration our future aim to apply the model to simulate the impact of the time evolution of the early biological pump (at the Neoproterozoic-Palaeozoic transition). For this reason, we don't push the range as far as what would 20 be needed to consider the impact of fast sinking rates typical of silicified or calcified small phytoplankton (McDonnell and Buesseler, 2010;Lam et al., 2011). In our sensitivity simulations, zrem S is in the range of [0, 40 m], while zrem L varies in the range of [50 m, 450 m].

Timescales 25
Starting from the initial values listed in Table 3, the modelled state variables evolve towards equilibrium for any couple of values of zrem S and zrem L in the explored interval. Simple mass conservation checks show no hidden source or sink of tracers in the model's boxes. Figure 3 illustrates an example of evolution of the variables for zrem S and zrem L in the middle of the interval of explored values for both particle types. In all the ocean boxes, P shows an initial oscillation that evolves on timescales of tens of thousands of years (Figure 3a,b), as expected by 30 the typical timescale of evolution of the tracer (Lenton and Watson, 2000). This is followed by a slower drift which depends on the dynamics of the deep water oxygen content, as the release and burial of P in the sediments depends on the level of oxygenation of the deep ocean and especially of the deep shelf sea. P reaches complete equilibrium as soon as the deep ocean boxes become stably oxygenated. The timescales of evolution of O are slower and lay on the order of tens of millions of years (Lenton and Watson, 2000). Oxygen in the deep shelf overcomes hypoxia after the first few millions of years and then slowly evolves towards equilibrium on the same timescale of O in the other ocean boxes. The dynamics of SedPorg is also strongly driven by level of oxygenation of the deep shelf sea. 5

Modern ocean budgets and fluxes
Modern estimates of the zrem S and zrem L vary depending on the region of sampling and on the local community structure, with most of the measurements focusing on large or heavy particles and most studies focusing on the open ocean (Iversen and Ploug, 2010;Cavan et al., 2017;Lam et al., 2011). Furthermore, only a very limited number of measurements accounts for both microbial and zooplankton remineralization, the latter disregarded by 10 lab measurements of zrem . Considering the fundamental role of the shelf sea in our model (always accounting for > 98 % of the total burial), we evaluate modelled tracer budgets and fluxes for values of zrem L that lay around 76 m, as measured in situ by Cavan et al. (2017) for a modern shelf sea. We pose no restrictions on zrem S due to the lack of precise measurements. A summary of our evaluation is provided in Table   5. 15 In the above mentioned range of zrem, our model predicts equilibrium budgets of between 3200 TmolP and 3400 TmolP for phosphorus, and an oxygen budget of between 100 PmolO2 and 150 PmolO2 in the entire ocean, compared to the estimated total P reservoir of 3100 TmolP (Watson et al., 2017) and estimated deep ocean O2 reservoir of 220 PmolO2 (Slomp and Van Cappellen, 2006). Due to the relative size of the ocean boxes, it is important to underline that total budgets are strongly driven by the deep open ocean budget. Modelled biological fluxes, such as production and export, fall around the low end of present estimates ( Figure   7). Our model predicts an equivalent C primary production of between 1400 TmolC yr -1 and 3000 TmolC yr -1 , and an export below the euphotic layer ranges between 300 TmolC yr -1 and 430 TmolC yr -1 , both calculated using a fixed C:P Redfield ratio of 106. These must be compared to an expected value of production of at least 3300 35 TmolC yr -1 (Carr et al., 2006) and an estimated export flux of at least 415 TmolC yr -1 (Henson et al., 2011). Despite the absolute fluxes being lower than expected, modelled export production (the export to production ratio) and the burial to production ratio compare well to range of present estimates. The modelled export corresponds to between 10 % and 32 % of total production, strongly depending on zrem S , compared to an expected range of 2 % - 20 % (Boyd and Trull, 2007). Buried Porg corresponds to between 0.3 % and 0.7 % of total production, compared to an expected 0.4 % (Sarmiento and Gruber, 2006).
In terms of the shelf contribution to the total fluxes, model results also fall close to present estimates. Modelled production in the surface shelf sea represents between 12% and 20% of total production (expected 20%) (Barrón and Duarte, 2015;Wollast, 1998). The fraction of modelled export and burial that happens in the shelf region 5 represent, respectively, [20 %, 23 %] and nearly 100 % of the total ocean fluxes, compared to estimated modern values of 29 % and 91 % (Sarmiento and Gruber, 2006)

Budgets and fluxes sensitivity to changes in zrem
Around the lowest values of zrem L adopted in the present study, i.e., in the range of [50 m, 100 m], our model shows a strong sensitivity of the total and local ocean P and O budgets for small changes of zrem L (Figure 4). This 15 is true for any zrem S , with minor differences between low and high zrem S values. For smaller zrem L , the model shows a sharp increase in P concentrations in all the ocean boxes and a substantial decrease of O levels at depth ( Figures   5,6), which are coupled to high levels of production and remineralization and low rates of sedimentation ( Figure   7). Essentially slow sinking and/or rapid remineralization results in inefficient removal of P to shelf sea sediments, requiring the ocean concentration of P to rise considerably for P output to balance (fixed) P input to the ocean. for nearly any value of zrem S and zrem L , the deep shelf is hypoxic or even suboxic for a broad range of small values 25 of zrem L , especially close to modern shelf zrem L observations. Considering the wide spatial extension of our boxes, we expect these low oxygen levels to indicate the development of local anoxia in the deep shelf.
In a limited interval of small zrem S values (roughly zrem S < 6 m), model results depend only on the LPorg properties due to the rather irrelevant contribution of SPorg to export and remineralization. For larger zrem values (zrem S > 6 m and zrem L > 100 m), model results show a strong interdependence of equilibrium budgets and absolute fluxes on 30 both zrem S and zrem L . Interestingly, in this range of values, export production depends nearly entirely on the small particle properties, ranging between 10 % for low zrem S and 32 % for high zrem S , an overall trend that affects also the ratio of deep remineralization to surface production (Figure 7).

Budgets and fluxes contribution by particle class
The relative role of small and large particles to modelled biological and physical fluxes depends on a combination of their inherent properties (zrem) and of coagulation. Coagulation of SPorg into LPorg after production in surface boxes affects between 13 % and 55 % of the small particles in the surface shelf sea and between 4 % and 28 % of the small particles in the surface open ocean (Figure 8a,b). The highest rates of coagulation of SPorg into LPorg in 5 the surface ocean are found for especially high P concentrations resulting in high rates of small particle production.
Rates decrease quickly moving away from these high P condition. These correspond to roughly zrem L < 100 m, meaning rather labile or light large particles, which contribute poorly to P removal and net O production.
For zrem L > 100 m, LPorg efficiently remove P from the water column, limiting production. Coagulation rates are therefore lower and vary in a more limited range of values. For this reason, the contribution of SPorg to the total 10 export below the euphotic layer is strongly dominated by the value of zrem S , with a null contribution to export for all values of zrem S < 10 m and increasing values above it. This trend is reflected in the deep-water small particle fraction (Figure 8c,d). Coagulation rates at depth, even though lower compared to the surface, still reach values of a few percent in the shelf sea (≤ 6 % of the exported SPorg). Small particles contribute up to 66 % to export and up to 25 % to the sediment accumulation in the shelf sea, with the highest contribution to sediment accumulation 15 being reached for large zrem S and low zrem L . In the open ocean, small particles represent up to 87 % of the total export, with the percentage being strongly dependent on the value of zrem S . Our model highlights therefore the different role of large and small particles in the determination of the equilibrium budgets and fluxes. Coagulation into large (fast sinking, less liable) particles is essential to maintain high enough sedimentation and burial rates, therefore allowing O accumulation in the system. At the same time, small (slow sinking, more liable) particles 20 tune the total magnitude of export and remineralization below the euphotic layer, impacting the distribution of oxygen and nutrients throughout the water column.

Box model limitations
BPOP consists in a simple box model with 4 ocean boxes, 2 sediment boxes and 1 atmospheric box. As with every box model, BPOP only allows a very rough and fundamental representation of the ocean's topography and circulation as well as of the exchange fluxes between ocean, atmosphere and sediments. Even though this may be a limitation in the context of the study of the well-known modern (and future) ocean, such a computationally 30 inexpensive model can be a useful tool to for a first exploration of a large variety of projected conditions. In the context of understanding past ocean changes, often characterized by a limited availability of observational data, the use of such a simple model constitutes instead an effective and honest approach to understand global shifts in budgets and fluxes. Furthermore, BPOP explicitly distinguishes between the well sampled shelf sea and the less known open ocean of deep time, therefore allowing to relate shelf data with large scale open sea conditions. 35

Limitations connected to the biological pump representation
In our model we adopt a very simplified representation of the biological pump, including two particle classes, "small" and "large", generated by production and coagulation, assuming that, on average, zrem S < zrem L . This scheme resembles the one commonly used in ocean biogeochemical models Jackson and Burd, 2015). Our model does not include a DOM pool for reasons mostly connected to the implicit representation 5 of the biological pump and the complete remineralization of the non-sedimented organic material at each integration step.
Modelled particles get remineralized through the water column according to their characteristic zrem. Even though for simplicity we do not use a continuum spectrum of zrem, the use of two particles classes is in line with observations showing two distinct peaks in the observed distribution of particles' sinking speeds (Riley et al., 10 2012;Alonso-González et al., 2010). Furthermore, this simplification still allows to closely approximate the empirical particle flux curve as a function of depth, also known as Martin's curve (Boyd and Trull, 2007).
We assume that zrem S and zrem L do not vary between the shelf sea and the open ocean. However, modern ocean observations show cross-shore changes in the phytoplankton community structure and sinking speeds (Barton et al., 2013). Our simplifying assumption may therefore cause the overestimation of the relative contribution of the 15 shelf sea to the total burial flux of Porg. Despite this, we believe that this choice is still convenient in the context of the current model, as it allows us to reduce the number of parameters in such a simple box model representation of the ocean's biological pump.
Observations suggest that hard shelled phytoplankton types, especially calcified cells, contribute substantially to the vertical export and burial of the organic material thanks to extremely large zrem despite their small size (Lam 20 et al., 2011;Iversen and Ploug, 2010). In the present study we focus on an interval of zrem S and zrem L values that are most likely to resemble the biological pump conditions of the Neoproterozoic -early Paleozoic ocean, before the evolution of such phytoplankton types. However, the model allows to explore different ranges of zrem S and zrem L values and to tune the rate of coagulation in order to explore the influence of these phytoplankton classes.
Even though bacterial remineralization is thought to be the dominant pathway for organic matter recycling on a 25 global scale, especially at low latitudes (Rivkin and Legendre, 2001), modern ocean coastal environments are also characterized by high grazing rates. The evolution of zooplankton and increasingly large grazers may have had a different impact on the effective zrem S and low zrem L , given additional Porg transformations such as particle fragmentation due to sloppy feeding Iversen and Poulsen, 2007). These processes can limit the large particle burial rates, while resulting in the deep production of small particle, s-POM and DOM. Our 30 model does not currently account for particle fragmentation, however the process could be easily considered in future model developments. In this context, new processes such as the sedimentation and burial of large grazers should also be considered.

Sensitivity to parameter choices
We discuss here the model sensitivity to changes in a set of significant parameters adopted to describe its 35 geometry, circulation and biological processes. Overall, none of the sensitivity experiments showed significant changes in the model results and conclusions: trends in budgets and fluxes obtained varying zrem S and zrem L , as well as our main results regarding the relative deep shelf and open ocean oxygenation remain unchanged.
Among the geometrical box model parameters, a key value is represented by the percentage of shelf sea area ( shelf). An increase (e.g., doubling) in shelf results in an overall decrease in the total budget of P and increase in O due to the larger ratio of burial to production, which is facilitated by a larger extension of the surface of shallow water. Interestingly, deep shelf anoxia is enhanced for larger shelf, i.e., anoxia is observed for a wider range of zrem S and zrem L values, while the deep ocean tends to be more oxygenated. Despite a doubling of shelf, however, 5 model results largely remain in the same range of those found for modern shelf.
We explored the effect of varying the physical circulation parameters. Changes in upwelling (Upw), have an important impact on the modelled ocean's budgets. An increase in Upw induces a lowering of P levels, especially in the deep shelf, due to their recirculation towards the surface and consequent uptake by production. This is Lastly, we explored the model sensitivity to the choice of key biogeochemical parameters representing rates of transformation. Both the use of higher coagulation rates (cgr) and the use of higher rates of formation of mineral 25 Ca-P (CaPr) result in a general increase in O levels and decrease in nutrient availability due to larger sedimentation and burial rates. However, we find again no substantial change in the model behaviour nor in the relative contribution to budgets and fluxes of each modelled ocean box.

Past changes in the biological pump 30
The evolution of larger and heavier cells during the Neoproterozoic and across the Neoproterozoic-Paleozoic transition is hypothesised to have caused significant changes in the ocean's nutrient and redox state ). Our new model can be used to assess the impact of this evolution in both the shelf and the open ocean. Our first model results highlight that for small zrem L , i.e., for an early biological pump with reduced capacity of export and burial, nutrient levels and production rates are particularly high. At the same time, an increase in 35 zrem S alone, fuelling higher remineralization rates at depth, can induce anoxia in the deep shelf while still maintaining the deep open ocean substantially oxygenated. The possibility of a coexistence of an anoxic deep shelf with an oxygenated deep open ocean has important implications for the interpretation of deep time redox proxy data, which come almost exclusively from shelf and slope environments, yet have been widely used to infer deep ocean anoxia for most of the Proterozoic Eon (Lenton and Daines, 2017) . We plan to use our model to further explore these changes in a time-frame perspective, introducing time varying boundary conditions (such as changes in Pin) and parameter properties.
Phytoplankton evolution as well as the development of heavier and larger marine organisms continued throughout the Phanerozoic (Katz et al., 2007). BPOP can also be used to explore the role of the biological pump in the onset 5 of OAEs in the course of the Mesozoic era, likely induced by enhanced productivity due to an upwelling intensification (Higgins et al., 2012). During the Mesozoic era, the evolution of dinoflagellates, calcareous and silica-encased phytoplankton also likely impacted the export and burial rates in a significant way (Katz et al., 2004). By extending the range of explored values of zrem S and zrem L , or possibly including the effect of grazing and/or an additional heavy POC class for shelled organisms, BPOP can also be used to study the consequences of 10 such evolution.

Future changes in the biological pump
Predicted future changes connected to global warming include, among the others, changes in ocean temperature, pH and stratification , with additional repercussions on plankton community structure, production, remineralization and export rates (Laufkötter et al., 2017;Acevedo-Trejos et al., 2014;Kwon et al., 15 2009). Our results show that around values of zrem L measured for modern shelf environment  modelled equilibrium budgets and fluxes are very sensitive to small changes in zrem. This indicates a potentially high sensitivity of the modern ocean to small changes in the biological pump, which may be particularly important in the deep shelf, where the boundary with suboxia is especially close (Keeling et al., 2010). Our model can be used to get a first assessment of the large-scale combined effect of predicted changes in the biological pump with 20 expected shifts in the physical ocean properties.

Exploring past and future changes in geometry, physics and biogeochemistry
In the present study we have focused on the impact of changes of zrem S and zrem L on the equilibrium budget and fluxes in the ocean. However, BPOP can be used to explore the effect of global changes in other physical or biogeochemical processes coupled to the biological pump dynamics. Aside from testing the robustness of our 25 results, the sensitivity tests presented in subsection 5.1.3 serve also as a first exploration of the possibility to apply the model to these further studies. We discuss here a few examples of past changes that could be explored with the present model. Through Earth's history, variations in the distribution of continents and in the mean sea level height likely impacted the percentage of shelf sea area ( shelf) throughout the global ocean (Katz et al., 2007). Changes in 30 climate and therefore in the mean temperature are expected to have affected both the air-sea gas exchange of oxygen (Schmidt number, ℎ ( )) and vertical mixing (Mixvo) (Petit et al., 1999). Reduced vertical mixing in warm periods is also expected to be relevant in the future because of global warming .
Changes in temperature are also known to impact biological activity directly, e.g,, by increasing remineralization rates (rmr) (Laufkötter et al., 2017), and indirectly, e.g., affecting production and mortality rates through changes 35 in the mixed layer depth (Polovina et al., 1995). Lastly, climatic shifts can also cause changes in the intensity of alongshore winds and therefore in the upwelling circulation (Sydeman et al., 2014).

Conclusions and Outlook
This paper provides a description, evaluation and discussion of the new BPOP model. BPOP is aimed at exploring the effects of changes in the biological pump on the shelf and open ocean nutrient and redox state as well as on P and O fluxes. This model can be adopted for a large variety of studies aimed at exploring the impact of changes in the biological pump, i.e., the particle remineralization length scale zrem, in past and future ocean settings.

Code availability
The code is available for download in the supplementary material of the present publication, which also includes the user's manual.

Author contributions 5
TL and EL conceived the study. EL conceived and implemented the model. EL and TL evaluated and improved the model. Both authors contributed to the interpretation of the results, and to the writing of the present manuscript.

Competing interests
The authors declare that they have no conflict of interest.       10 unmeasuredgiven the analogous adopted functional form, we assume Ca-P formation to happen on a timescale close to that of Porg coagulation in the water column.  Sarmiento and Gruber (2006) Shelf sea production 12 % -20 % 20 % of total Prod Barrón and Duarte (2015); Wollast (1998) Shelf sea export 20 % -23 % 29 % of total Export Sarmiento and Gruber (2006) Shelf sea burial 100 % 91 % of total Burial Sarmiento and Gruber (2006)      (b) Porg export through the euphotic layer depth; (c) Export production, i.e. export to production ratio (d) Burial to production ratio.