Articles | Volume 13, issue 4
Model description paper
09 Apr 2020
Model description paper |  | 09 Apr 2020

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

Elisa Lovecchio and Timothy M. Lenton

The biological pump of the ocean has changed over Earth's 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 hypothesized to have caused important changes in the nutrient and redox properties of the ocean. To explore these hypotheses, we present here a new box model including oxygen (O), phosphorus (P) and a dynamical 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. Modeled O and P budgets and fluxes lie reasonably 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 the inefficient removal of P 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 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 with the biological pump in a variety of case studies.

1 Introduction

The “biological pump” describes the production of organic matter at the surface of the ocean (an oxygen source), its downward export/sinking flux, remineralization at depth (an oxygen sink) and burial. This set of processes acts against the homogenization of tracer concentrations by the ocean 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 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<60mmol m−3), suboxic (O2<5mmol m−3) or even anoxic – as 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 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 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).

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 1000 to 541 million years ago, as being 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 oxygen in the atmosphere (of uncertain cause) drove the oxygenation of the deep ocean at this time through air–sea 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 inventory (Lenton and Daines, 2018).

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 hypothesized 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 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 modeling work has examined the impact of changes in the organic matter remineralization length/depth (zrem) in the 3-D 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 concentrations – as it gets larger the oxygen 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 3-D 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 modulating 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 idealized 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 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. The model has a deliberately simplified treatment of redox carriers and is designed to focus on ocean P and ocean redox steady states, not on longer-term controls on atmospheric oxygen. In the following sections we describe the model, provide an evaluation of its performance in the context of modern observations and flux estimates, and finally present and discuss our model results.

Table 1List of the state variables in the model and of their units.

Download Print Version | Download XLSX

Table 2Parameters set that describes geometry of the box model.

Notes: (1) we assume a constant average euphotic layer depth of 100 m in both shelf and open sea; (2) the shelf sea is assumed to be 200 m deep in total, in line with the definition of shelf sea by Barrón and Duarte (2015); and (3) we assume an average open ocean depth of 3600 m (including euphotic layer).

Download Print Version | Download XLSX

2 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, 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 manual (see the Supplement) for further information on how to run the model.

2.1 Variables and circulation

The box model resolves explicitly for each relevant box the local concentrations of three types of tracers: molecular oxygen O2 (O), inorganic dissolved phosphorus (P) and sediment organic phosphorus (SedPorg). The total budgets of P and O, respectively PTOT and OTOT, 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 state and diagnostic variables in the model and their units are listed in Table 1. In the following subsections we describe the geometry of the box model and discuss the physical and geochemical fluxes that drive the dynamics of the tracers. Box properties are listed in Table 2, while the set of parameters adopted for the modeled physical and geochemical fluxes can be found in Table 3.

Table 3Parameter set pertaining to the initial conditions of the model, circulation mass fluxes and boundary fluxes.

Notes: (1) Chavez and Messié (2009) estimate 5.5 Sv (sverdrups) in the four major upwelling systems alone; (2) compare to 38 Sv (Sarmiento and Gruber, 2006), 17 Sv of mixing flux in the Southern Ocean alone (Meyer et al., 2015), estimated open ocean downwelling 38.5 Sv and upwelling 34.5 Sv (Ganachaud and Wunsch, 2000); (3) cross-shelf mass exchange due to lateral recirculation, tides and mixing aimed at including exchange processes other than upwelling (Fennel et al., 2005; Cole et al., 2015; Wollast, 1998); (4) minimal assumption for vertical mixing in nearshore regions due to seasonal and eddy mixing; see also Sect. 3.2 Sensitivity to parameter choices; (5) up to 70 % of river outflow reaches the open ocean; see Sharples et al. (2017); and (6) calculated from the equilibrium solution given Pin.

Download Print Version | Download XLSX

2.1.1 Box properties and physical fluxes of inorganic tracers

The box model includes four ocean boxes, one atmospheric box and two sediment boxes (Fig. 1a). The ocean and sediment boxes are equally split between shelf sea and open ocean, both including one surface ocean box and one deep ocean box.

O and P are exchanged between the four 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 ith box, the physical exchange flux (in millimoles per cubic meter per year) is represented by

(1) AdvMix ( C ) i = j MassFlux i j / V i ( C j - C i ) ,

where MassFluxij represents the volumetric flow between the ith box and any adjacent box j, while Vi is the volume of the ith box.

Figure 1Box model scheme with a representation of the physical and boundary fluxes affecting inorganic tracers in the water column and atmosphere, where blue arrows indicate advective and mixing fluxes and yellow arrows indicate air–sea gas exchange fluxes. The model includes seven boxes: surface shelf (ss), deep shelf (ds), surface open ocean (so), deep open ocean (do), atmosphere (at), shelf sediments (s) and open ocean sediments (o).


Figure 2Representation of the physical and biogeochemical fluxes affecting the Porg cycling in the model. Even though some processes (such as burial as Ca–P) are here represented in detail only in one box, the set of biogeochemical processes regulating the Porg dynamics in shelf sea and open ocean (both water column and sediments) is the same, as described in Sect. 2.2.


For each surface box i, air–sea gas exchange allows O fluxes between the ocean and the atmosphere (at). The flux (in millimoles of O2 per cubic meter per year) is positive when directed into the ocean and depends on the gas transfer velocity KW, atmospheric pressure pat (here assumed constant) and Henry's constant KHenry, as in

(2) AirSea i = K W ( O at p at / K Henry - O i ) A i / V i ,

where KW (in meter per year) is a function of the prescribed mean temperature Tmean and wind speed Wspeed (Sarmiento and Gruber, 2006).

2.1.2 Initialization and boundary fluxes

The model is initialized with an even concentration of P (Pini) in all the ocean boxes, zero oxygen (Oini) and zero SedPorg ((Porg)ini). 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 organic phosphorus (Porg) burial flux balances Pin. Oxidative weathering determined by atmospheric oxygen Oat constitutes a net sink flux for O. The weathering flux (per year) depends on a constant baseline flux W0, and it scales like the square root of the oxygen mixing ratio normalized to present values Omix0 (Lenton et al., 2018), following

(3) OxyWeath = W 0 O at / O mix 0 .

2.2 Biological pump details

The modeled tracer cycles are coupled with a set of biological transformations, i.e. the biological pump, governing the cycle of production, and 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 water column, and we only calculate fluxes of water column Porg without treating Porg as a state variable. This scheme is similar to the one used to represent detrital POM in some modern ocean biogeochemical models (Moore et al., 2004). P and O biological fluxes are coupled with 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.

Table 4Parameters set pertaining to the Porg cycle of the model and coupled biogeochemical fluxes.

Notes: (1) maximum P uptake rate, meant to account for environmental limitations of phytoplankton growth rate (such as light and temperature); the magnitude of the rate takes into account that we are not explicitly resolving phytoplankton concentrations (order of 10−2mmol P); see also production in Gruber et al. (2006) and Yool and Tyrrell (2003); (2) measured values vary in the range from 0.01 mmol m−3 up to a few millimoles per cubic meter, varying for different phytoplankton types; see Lomas et al. (2014), Tantanasarit et al. (2013), Krumhardt et al. (2013), Lin et al. (2016) and Klausmeier et al. (2004); (3) measured half-saturation constant for oxygen uptake varies in the range of 0.1–3 mmol m−3 (Ploug, 2001); (4) biogeochemical models commonly switch to anaerobic respiration below 4 mmol m−3 (Paulmier et al., 2009); measurements suggest a value close to 19 mmol m−3 (DeVries and Weber, 2017); (5) Cavan et al. (2017) showed that small particles are about 85 % of the total sinking particles abundance in the coastal region at export depth; the parameter was further tuned to bring the model closer to modern ocean conditions; (6) on the same order of magnitude as Gruber et al. (2006); (7) unmeasured – given 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 in models with explicit particle pools (Gruber et al., 2006).

Download Print Version | Download XLSX

2.2.1 Particle classes, production and coagulation

The model includes two Porg classes that get produced, exported and remineralized in the ocean water column, which are 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 a better reproduction of 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).

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 millimoles of P per cubic meter per year) in each ith box only generates SPorg, according to

(4) Prod i = P eff ( P i / ( P i + K P ) ) P i .

LPorg is generated via the coagulation of SPorg at the surface after production. As we do not explicitly solve for the concentrations of SPorg and LPorg, we assume that the coagulation (in millimoles of P per cubic meter per year) of SPorg into LPorg in each box i is proportional to the rate of production of small particles as follows:

(5) Coag i = cg f Prod i .

This is a necessary simplifying assumption compared to the usual coagulation models, which define the flux as the square of the particle concentration (Boyd and Trull, 2007; Gruber et al., 2006), given the fact that our model does not resolve this variable. 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.

2.2.2 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 paper, SPorg and LPorg are redistributed throughout the water column exclusively by implicitly modeled gravitational sinking before being buried, accumulated in the sediments or remineralized. Even though the vertical export by 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 (zremS=0), these fluxes are not currently accounted for in the model.

2.2.3 Remineralization length scheme

The export and sedimentation fluxes of Porg through the water column are represented by a remineralization length 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 as follows:

(6) f k ( z ) = f 0 k e - z - z 0 z rem k ,

where f0k is the flux (in millimoles of P per square meter per year) at the reference depth z0, and the index k 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 choice of z0 (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 meters and can be calculated as the ratio between the particle sinking speed and the remineralization rate of the particle (Cavan et al., 2017). Consequently, zrem implicitly contains information on several particle-inherent properties (among which are density, size, shape and organic matter liability) as well as information about the surrounding environment, e.g. the type of heterotrophs that 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.

2.2.4 Sediments and burial

SPorg and LPorg accumulate in the sediments as SedPorg, which is calculated as a density per unit of area. The flux (in millimoles of P per square meter per year) into the sediment box i depends on the organic matter fluxes into 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 modeled as a function of the square of SedPorg that accumulates in the sediments and is regulated by a constant rate coefficient CaPr. Ca–P formation happens at a lower rate under low-oxygen conditions (CaPr=CaPrfsan with fsan<1), in agreement with observations and previous models (Slomp and Van Cappellen, 2006). The transition 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 (in millimoles per square meter per year) as in


This flux is essential to balance the continuous P river input, therefore preventing the ocean from overflowing with nutrients.

2.2.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. The overall loss of P due to export (in millimoles of P per cubic meter per year) from a surface box i to a deep box j via gravitational sinking, is calculated as

(9) VExp i = ( Flx_SP org i exp ( - ( Δ Z eu / 2 ) / z rem S ) + Flx_ LP org i exp ( - ( Δ Z eu / 2 ) / z rem L ) ) / Δ Z i ,

where the fluxes per unit of area of SPorg and LPorg in the surface boxes depend on production and coagulation as described in Sect. 2.2.1.

At depth, the remineralization of Porg that does not reach the sediments happens through both aerobic and anaerobic processes, completely depleting the remaining Porg. The amount of inorganic P released in each deep box j by water-column remineralization (in millimoles of P per cubic meter per year) is therefore calculated as


In each deep ocean box i, aerobic remineralization uses some of the available oxygen and is therefore limited by Michaelis–Menten kinetics with a half-saturation constant KwO (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. In our model, the reducing agent produced by anaerobic remineralization is methane gas, and it is only produced when the sediments and the deep shelf water column have gone anoxic. As we do not track other oxidizing agents such as SO4, there is nothing for the methane to be oxidized by until it reaches the surface ocean, and as the surface ocean is equilibrated with the atmosphere, the fact that we assume oxidation in the atmosphere is a reasonable approximation. 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 rather 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 atmosphere and happens at a faster rate rmr=rmrfean with fean>1, in agreement with recent observations and previous models (Slomp and Van Cappellen, 2006). The release of Porg from a sediment box i into the overlaying ocean box due to sediment remineralization (in millimoles of P per cubic meter per year) is therefore the sum of the two terms as in


2.3 Equations summary

The dynamics of the 11 state variables in the model are 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 dimensions of the reference boxes or number of moles (atmosphere). A full set of equations including the explicit formulation of all the flux terms for each box can be found in the Appendix.


where Pin is the river input of P to the surface of the ocean and AdvMix indicates the advective and mixing physical fluxes of the variable of interest (which differ for each box according to the circulation scheme); Exp is the export flux of Porg in P units; WcRem indicates the water column complete remineralization of the organic material in P units, which is split into an anaerobic (Ana) and aerobic (Aer) component; SedRem indicates the sediment remineralization of SedPorg in P units (also aerobic and anaerobic); AirSea represents the air–sea flux exchange of O; OxyWeath is the O weathering flux sink; SedFlx is the SedPorg accumulation flux as regulated by the remineralization length scheme at the bottom of the water column; and finally CaPform represents the sediment burial flux of P in mineral form. For each box, flux terms are rescaled with the appropriate box geometry.

2.4 Strategy: sensitivity studies for varying zrem

In order to characterize the model, we analyze the equilibrium budgets and fluxes of the state variables for varying zrem values separately for SPorg and LPorg, respectively zremS and zremL. We adopt a range of zrem values that fall close to modern observations (Cavan et al., 2017; Buesseler and Boyd, 2009; Marsay et al., 2015) and takes 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–Paleozoic transition). For this reason, we do not push the range as far as what would 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, zremS is in the range of 0–40 m, while zremL varies in the range of 50–450 m.

3 Evaluation

3.1 Timescales

Starting from the initial values listed in Table 3, the modeled state variables evolve towards equilibrium for any pair of values of zremS and zremL in the explored interval. Simple mass conservation checks show no hidden source or sink of tracers in the boxes of the model. Figure 3 illustrates an example of evolution of the variables for zremS and zremL 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 (Fig. 3a, b), as expected by 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 depend 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 the evolution of O are slower and lie 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 are also strongly driven by level of oxygenation of the deep shelf sea. The dynamical response of the model to changes in the biological pump is rapid subsequent to the model equilibrating considering the given initial conditions. For example, step changes in the zrem values of the particles result in a transition time to a new equilibrium that is on the order of a few tens of thousands of years, which is the typical timescale of the P cycle.

Figure 3Evolution of the state variables from the initial conditions listed in Table 2 and remineralization lengths roughly in the middle of the interval of explored values, zremS=20m and zremL=250m. (a) Evolution of inorganic phosphorus P in the water column (left axis) and of organic phosphorus in the sediments SedPorg (right axis). (b) A zoomed-in view of the dynamics of P in the first 200 000 years. (c) Evolution of oxygen in the water column (left axis) and atmosphere (right axis). In (c) the two lines Oss and Oso are overlapping; the two variables evolve closely due to the coupling of the surface ocean with the atmosphere via air–sea gas exchange.


3.2 Modern ocean budgets and fluxes

Modern estimates of the zremS and zremL 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 account for both microbial and zooplankton remineralization, the latter disregarded by lab measurements of zrem (Cavan et al., 2017). Considering the fundamental role of the shelf sea in our model (always accounting for >98 % of the total burial), we evaluate modeled tracer budgets and fluxes for values of zremL that lie around 76 m, as measured in situ by Cavan et al. (2017) for a modern shelf sea. We pose no restrictions on zremS due to the lack of precise measurements. A summary of our evaluation is provided in Table 5.

Table 5Summary of the model evaluation provided in Sect. 3. Modern observations and estimates are compared to model results obtained for zremL in the range of measured values for a modern shelf sea (Cavan et al., 2017).

Download Print Version | Download XLSX

In the above-mentioned range of zrem, our model predicts equilibrium budgets of between 2250 and 2970 Tmol P for phosphorus and an oxygen budget of between 100 and 107 Pmol O2 in the entire ocean, compared to the estimated total P reservoir of 3100 Tmol P (Watson et al., 2017) and estimated ocean O2 reservoir of between 225 and 310 Pmol O2 (Keeling et al., 1993; Duursma and Boisson, 1994). 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 and that the low-oxygen reservoir of our model may be connected to an underestimation of the deep open ocean oxygenation.

Deep shelf P and O concentrations lie in the ranges of 3.9–4.9 mmol m−3 and 3.8–9.2 mmol m−3, respectively (Figs. 5 and 6). Deep shelf nutrient concentrations are higher than expected by about a factor of 2 compared to modern values, possibly due to the fact that our model does not store any Porg in the water column or due to an underestimation of the vertical supply of nutrients to the surface shelf (e.g. via mixing). Limiting deep P concentrations via lower remineralization or higher burial rates, however, also results in sensibly lower production rates. In the deep open ocean, P and O concentrations fall in the ranges of 1.9–2.5 mmol m−3 and 76–83 mmol m−3, respectively. For any combination of zremS and zremL, O levels in surface ocean boxes lie between 273 and 274 mmol m−3, a good approximation of average modern surface values (Garcia et al., 2018b). In general, the deep shelf always shows the highest P values and lowest O concentrations compared to the other ocean regions, while, as expected, the surface shelf sea is richer in P compared to the surface open ocean.

Figure 4Total ocean budgets of (a) P and (b) O at equilibrium for varying zremS and zremL.


In order to compare the modeled fluxes to modern estimates, we converted our results into carbon (C) units assuming a C:P Redfield ratio of 106. However, recent studies found a substantially higher mean C:P ratio for the modern ocean (Martiny et al., 2014); therefore our derived C fluxes may be a conservative estimate. Modeled biological fluxes in C units, such as production and export, fall just below the low end of present estimates (Fig. 7). Our model predicts a total primary production of between 11 and 30 Gt C yr−1, and an export below the euphotic layer ranges between 3.4 and 3.8 Gt C yr−1. These must be compared to an expected value of production of between 35 and 80 Gt C yr−1 (Carr et al., 2006) and an estimated export flux of at least 4 Tmol C yr−1 (Henson et al., 2011). Despite the absolute fluxes being at the low end of the present estimates, modeled export production (the export to production ratio) and the burial to production ratio compare well to range of present estimates. The modeled export corresponds to between 11 % and 33 % of total production, strongly depending on zremS, compared to an expected range of 2 %–20 % (Boyd and Trull, 2007). Buried Porg corresponds to between 0.3 % and 1 % 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. Modeled production in the surface shelf sea represents between 16 % and 27 % of total production (expected 20 %) (Barrón and Duarte, 2015; Wollast, 1998). The fraction of modeled export and burial that happens in the shelf region represent, respectively, [16 %, 27 %] and nearly 100 % of the total ocean fluxes, compared to estimated modern values of 29 % and 91 % (Sarmiento and Gruber, 2006). Our overestimation of the shelf contribution to the burial fluxes may be due to the underestimation of the zrem of open-ocean particles compared to observations (Cavan et al., 2017; Lam et al., 2011), i.e. our choice of using the same value of zremS and zremL for both the coastal and the open ocean box. This simplifying assumption limits the capacity of Porg to reach the deep sediment layer in the open ocean. We explore potential limitations of this choice in the Discussion section.

4 Results

4.1 Budgets and fluxes sensitivity to changes in zrem

Around the lowest values of zremL adopted in the present study, i.e. in the range of 50–100 m, our model shows a strong sensitivity of the total and local ocean P and O budgets to small changes in zremL (Fig. 4). This is true for any zremS, with minor differences between low and high zremS values. For smaller zremL, the model shows a sharp increase in P concentrations in all the ocean boxes and a substantial decrease in O levels at depth (Figs. 5 and 6), which are coupled with high levels of production and remineralization and low rates of sedimentation (Fig. 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.

Figure 5Local P concentration in each ocean box for varying zremS and zremL for the following: (a) surface shelf sea, ss; (b) surface open ocean, so; (c) deep shelf sea, ds; and (d) deep open ocean, do. Surface ocean boxes, as well as deep ocean boxes, are plotted on the same scale.


Figure 6O concentrations at equilibrium for varying zremS and zremL for the following: (a) deep shelf sea, ds; (b) deep open ocean, do. Surface ocean boxes (not shown) have nearly constant values of O for any set of zrem due to the air–sea gas exchange, which strongly couples them to the atmosphere.


Our model results show that for any pair of values of zremS and zremL in the entire explored range, the biological pump is able to oxygenate the surface ocean (surface O levels lie close to 273 mmol m−3) and, for most values, also to maintain the deep ocean above the level of hypoxia (Fig. 6). The model shows a substantial difference between the deep shelf and the deep open ocean; while the latter is substantially oxygenated (O>50mmol m−3) for nearly any values of zremS and zremL, the deep shelf is hypoxic or even suboxic for a broad range of small values of zremL, especially close to modern shelf zremL 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 zremS values (roughly zremS<6m), model results depend only on the LPorg properties due to the rather irrelevant contribution of SPorg to export and remineralization. For larger zrem values (zremS>6m and zremL>100m), model results show a strong interdependence of equilibrium budgets and absolute fluxes on both zremS and zremL. Interestingly, in this range of values, export production depends very strongly on the small particle properties, ranging between 10 % for low zremS and 30 % for high zremS, an overall trend that also affects the ratio of deep remineralization to surface production (Fig. 7).

Figure 7Biological pump fluxes in P units for the entire ocean for varying zremS and zremL in the following cases: (a) Porg surface production; (b) Porg export through the euphotic layer depth; (c) export / production, i.e. export to production ratio; and (d) burial to production ratio.


It is also important to notice that, for any couple of zremS and zremL, modeled tracer concentrations and fluxes fall in a range of values that never exceeds by orders of magnitude the modern observed values. Considering all of the ocean boxes, P concentrations vary in the range of roughly 0.2 and 9 mmol m−3, while O levels lie between 0.5 and 205 mmol m−3. Production in carbon units lies in the interval [7.6, 70.7 Gt C yr−1].

4.2 Budgets and fluxes contribution by particle class

The relative role of small and large particles in modeled biological and physical fluxes depends on a combination of their inherent properties (zrem) and coagulation. In our simple model, coagulation of SPorg into LPorg after production in surface boxes affects a constant fraction (cgf=0.22) of the produced particles. This fraction was determined by model tuning to modern ocean conditions and lies close to modern ocean observations of the large-particle fraction (15% of the total particles) at export depth (Cavan et al., 2017).

For zremL>100m, LPorg efficiently removes P from the water column, limiting production. The contribution of SPorg to the total export below the euphotic layer, however, is strongly dominated by the value of zremS, with a null contribution to export for all values of zremS<10m and increasing values above it. This trend is reflected in the deep water small-particle fraction (Fig. 8c, d). Small particles contribute up to 73 % to export in both ocean boxes and up to 60 % to the sediment accumulation in the shelf sea, with the highest contribution to sediment accumulation being reached for large zremS and low zremL. 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 and 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 and more liable) particles tune the total magnitude of export and remineralization below the euphotic layer, affecting the distribution of oxygen and nutrients throughout the water column.

5 Discussion

5.1 Model limitations and robustness

5.1.1 General limitations

BPOP consists of a simple box model with four ocean boxes, two sediment boxes and one atmospheric box. As with every box model, BPOP only allows a very rough and fundamental representation of the topography and circulation of the ocean 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 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 understanding 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 the relation of shelf data with large-scale open ocean conditions.

The model deliberately simplifies the redox carriers and processes represented, neglecting denitrification and iron and sulfate reduction. Including additional oxidants and/or methane consumption in deeper water column would be expected to intensify anoxia results at depth. However, our current results suggest that the model is underestimating the ocean total oxygen budget overall, mostly driven by the deep open ocean reservoir. This suggests that neglecting these additional processes in our simple box model does not lead to an overestimation of oxygen accumulation at depth. Including additional state variables and processes could also lead to more complex dynamical behaviors (Wallmann, 2010).

We include anaerobic remineralization of Porg as being faster than aerobic degradation, but in reality this is not the case for carbon – which is remineralized at a similar or slower rate under anoxic vs. oxic conditions (Burdige, 2007; Hedges et al., 1999; Dale et al., 2015). Hence, in reality, under anoxic conditions, there is preferential regeneration of phosphorus and organic C:P burial ratios rise considerably, altering the long-term steady state of atmospheric oxygen (Van Cappellen and Ingall, 1996). We do not consider these aspects here because to do so would require adding state variables for organic carbon (as distinct from organic phosphorus) and because our focus here is on changes in ocean phosphorus and ocean redox under an unchanged oxygen steady state. In future work we intend to elaborate the model in order to explore long-term effects on atmospheric oxygen.

5.1.2 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, zremS<zremL. This scheme resembles the one commonly used in ocean biogeochemical models (Gruber et al., 2006; Jackson and Burd, 2015). Our model does not include a dissolved organic matter (DOM) pool for reasons mostly connected to the implicit representation of the biological pump and the complete remineralization of the non-sedimented organic material at each integration step. For the same reason, we do not resolve particle Porg concentrations, and therefore we model the coagulation flux as a constant fraction of production. A more physical representation of coagulation would require this flux to scale with the square of the particle concentrations (Boyd and Trull, 2007). Such a further development could potentially lead to increase large particle export for high surface P concentrations leading to high production and particle concentrations (and vice versa). We reserve this improvement as our first step for further model developments, which will include an explicit Porg representation.

Modeled particles get remineralized through the water column according to their characteristic zrem. Even though for simplicity we do not use a continuous spectrum of zrem, the use of two particle classes is in line with observations showing two distinct peaks in the observed distribution of the sinking speeds of the particles (Riley et al., 2012; Alonso-González et al., 2010). Furthermore, this simplification still allows us 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 zremS and zremL 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 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 biological pump of the ocean.

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 et al., 2011; Iversen and Ploug, 2010). In the present study we focus on an interval of zremS and zremL 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 the exploration of different ranges of zremS and zremL values and the tuning of 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 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 zremS and low zremL, given additional Porg transformations such as particle fragmentation due to sloppy feeding (Cavan et al., 2017; Iversen and Poulsen, 2007). These processes can limit the large-particle burial rates, while resulting in the deep production of small particles, suspended particulate organic matter (s-POM) and DOM. Our 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.

5.1.3 Sensitivity to parameter choices

We discuss here the model sensitivity to changes in a set of significant parameters adopted to describe its 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 zremS and zremL and 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 zremS and zremL values, while the deep ocean tends to be more oxygenated. Despite a doubling of 𝒫shelf, however, 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 budgets of the modeled ocean. 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 coupled with an overall larger equilibrium O budget due to higher storage in the deep open ocean and consequent recirculation into the deep shelf. Deep shelf suboxia is still possible but for a more limited range of zremL values. Changes in vertical mixing in the open ocean (Mixvo) affect the overall P and O budgets mostly for high zremL. For lower Mixvo, the O budget decreases due to lower O storage at depth, while P increases. Changes in vertical mixing on the shelf (Mixvs), instead, have a minor impact on the total budgets and fluxes of the model, while locally modulating shelf oxygen and nutrient concentrations. Lateral mixing fluxes (Mixls, Mixld) were included in our model for means of generalization and in order to account for the influence of non-upwelling margins, with a lower value than in previous studies (Fennel et al., 2005). Changes in Mixls and Mixld result in significant changes in the deep ocean storage of tracers and open ocean production, with little impact on the budget of the other ocean boxes. However, in this case also, our main conclusions remain unaffected.

We explored the impact of changing the portion of nutrients delivered directly to the open ocean, Popen. Even large changes in this parameter do not significantly affect the results of the model, indicating that the relative levels of P and O at equilibrium are determined by the internal physical and biogeochemical dynamics of the model, rather than by boundary conditions.

Lastly, we explored the model sensitivity to the choice of key biogeochemical parameters representing rates of transformation. Both increasing coagulation (cgf) and the use of higher rates of formation of mineral 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 behavior nor in the relative contribution to budgets and fluxes of each modeled ocean box.

Furthermore, we have tested the impact of having sediment remineralization rates that vary with the zrem of the particles, under the assumption that the liability of small and large particles may be different. In our experiment, we increased the remineralization rmr rate linearly with zrem by 40 % of our baseline value (rmr0), with rmr0 being found at the center of the interval of explored values of zrem=[0, 450 m]. Under these conditions, we obtained a higher decoupling between the influence of zremS and zremL on budgets and fluxes, both being more strongly driven by the small particle properties for large values of zremL.

5.2 Model applications

5.2.1 Past changes in the biological pump

The evolution of larger and heavier cells during the Neoproterozoic and across the Neoproterozoic–Paleozoic transition is hypothesized to have caused significant changes in the nutrient and redox state of the ocean (Lenton and Daines, 2018). 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 zremL, 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 zremS alone, fueling 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 from 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 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 and 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 zremS and zremL or possibly including the effect of grazing and/or an additional heavy POM class for shelled organisms, BPOP can also be used to study the consequences of such an evolution.

5.2.2 Future changes in the biological pump

Predicted future changes connected to global warming include, among others, changes in ocean temperature, pH and stratification (Gruber et al., 2004), 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., 2009). Our results show that around values of zremL measured for the modern shelf environment (Cavan et al., 2017), modeled 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 expected shifts in the physical ocean properties.

5.2.3 Exploring past and future changes in geometry, physics and biogeochemistry

In the present study we have focused on the impact of changes in zremS and zremL 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 with the biological pump dynamics. Aside from testing the robustness of our results, the sensitivity tests presented in Sect. 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 climate and therefore in the mean temperature are expected to have affected both the air–sea gas exchange of oxygen – Schmidt number NSch(Tmean) – 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 (Gruber et al., 2004). 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 in the mixed layer depth (Polovina et al., 1995). Climatic shifts can also cause changes in the intensity of alongshore winds and therefore in the upwelling circulation (Sydeman et al., 2014). Lastly, the model can be used to test the impact of changes in the biogeochemical cycles, including shifts in the Redfield ratio as well as global changes in the P input (Pin) to the ocean (Reinhard et al., 2017; Filippelli, 2008).

6 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. Furthermore, it allows us to couple changes in POM properties with changes in the geometry, circulation and boundary conditions of the ocean.

Despite its simple representation of the ocean circulation and the biological pump, the model can reasonably simulate values of the current P and O tracer budgets and biological pump fluxes. The model predicts potentially large variations in these P and O budgets and fluxes for past and future changes in the POM remineralization length. Our preliminary results also indicate that the early ocean may have been nutrient rich, with high levels of production and remineralization and that a suboxic deep shelf setting may have been compatible with an oxygenated deep open ocean.

We plan to apply this model to study the time evolution of the P and O budgets in both the shelf and the open ocean environment across the Neoproterozoic–Phanerozoic transition. Further developments of the model will be aimed at accounting for successive evolutionary innovations, including particle fragmentation due to grazing.

Appendix A: Equations

A1 Air–sea gas exchange of oxygen


A2 Surface shelf sea (ss)


A3 Deep shelf sea (ds)


A4 Surface open ocean (so)


A5 Deep open ocean (do)

(A36)Vdo=ΔZdoAocean(1-Pshelf)(A37)VIn_SPorgdo=VExp_SPorgso(A38)VIn_LPorgdo=VExp_LPorgsoRem_SPorgdo=VIn_SPorgdo(A39)1-exp-ΔZdo/zremS/ΔZdoRem_LPorgdo=VIn_ LPorgdo(A40)1-exp-ΔZdo/zremL/ΔZdoAerRem_SedPorgdo=rmrSedPorgo/ΔZdo(A41)Odo/Odo+KOsAnaRem_SedPorgdo=(rmrfean)SedPorgo/ΔZdo(A42)1-Odo/Odo+KOsdPdodt=(UpwPso-Pdo+MixldPds-Pdo+MixvoPso-Pdo)spy/Vdo+(Rem_SPorgdo+Rem_LPorgdo+AerRem_SedPorgdo(A43)+AnaRem_SedPorgdo)AerRemWcOdo=OPRedRem_SPorgdo+Rem_LPorgdo(A44)Odo/Odo+KOw(A45)AerRemSedOdo=OPRedAerRem_SedPorgdodOdodt=(UpwOso-Odo+MixldOds-Odo+MixvoOso-Odo)spy/Vdo-AerRemWcOdo(A46)-AerRemSedOdo

A6 Shelf sea sediments (s)


A7 Open ocean sediments (o)


A8 Atmosphere (at)


A9 Diagnostics: total budgets of P and O

Code availability

The code is available for download in the Supplement of the present publication, which also includes the user manual.


The supplement related to this article is available online at:

Author contributions

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 paper.

Competing interests

The authors declare that they have no conflict of interest.


We would like to thank Richard Boyle for his valuable input and suggestions. We further acknowledge the precious input of the three referees and the editor, which substantially improved the model and paper. This research was funded by NERC in the framework of the project Biosphere Evolution, Transitions and Resilience (BETR).

Financial support

This research has been supported by the Natural Environment Research Council (grant no. NE/P013651/1).

Review statement

This paper was edited by Andrew Yool and reviewed by Klaus Wallmann and one anonymous referee.


Acevedo-Trejos, E., Brandt, G., Steinacher, M., and Merico, A.: A glimpse into the future composition of marine phytoplankton communities, Front. Marine Sci., 1, 15 pp.,, 2014. 

Alonso-González, I. J., Arístegui, J., Lee, C., Sanchez-Vidal, A., Calafat, A., Fabrés, J., Sangrá, P., Masqué, P., Hernández-Guerra, A., and Benítez-Barrios, V.: Role of slowly settling particles in the ocean carbon cycle, Geophys. Res. Lett., 37, L13608,, 2010. 

Anderson, L. A. and Sarmiento, J. L.: Redfield ratios of remineralization determined by nutrient data analysis, Global Biogeochem. Cy., 8, 65–80,, 1994. 

Baker, C. A., Henson, S. A., Cavan, E. L., Giering, S. L., Yool, A., Gehlen, M., Belcher, A., Riley, J. S., Smith, H. E., and Sanders, R.: Slow-sinking particulate organic carbon in the Atlantic Ocean: Magnitude, flux, and potential controls, Global Biogeochem. Cy., 31, 1051–1065, 2017. 

Barrón, C. and Duarte, C. M.: Dissolved organic carbon pools and export from the coastal ocean, Global Biogeochem. Cy., 29, 1725–1738,, 2015. 

Barton, A. D., Pershing, A. J., Litchman, E., Record, N. R., Edwards, K. F., Finkel, Z. V., Kiørboe, T., and Ward, B. A.: The biogeography of marine plankton traits, Ecol. Lett., 16, 522–534,, 2013. 

Boyd, P. and Trull, T.: Understanding the export of biogenic particles in oceanic waters: is there consensus?, Prog. Oceanogr., 72, 276–312, 2007. 

Boyle, R., Dahl, T., Bjerrum, C., and Canfield, D.: Bioturbation and directionality in Earth's carbon isotope record across the Neoproterozoic–Cambrian transition, Geobiology, 16, 252–278, 2018. 

Brocks, J. J., Jarrett, A. J., Sirantoine, E., Hallmann, C., Hoshino, Y., and Liyanage, T.: The rise of algae in Cryogenian oceans and the emergence of animals, Nature, 548, 578–581,, 2017. 

Buesseler, K. O. and Boyd, P. W.: Shedding light on processes that control particle export and flux attenuation in the twilight zone of the open ocean, Limnol. Oceanogr., 54, 1210–1232, 2009. 

Burdige, D. J.: Preservation of organic matter in marine sediments: controls, mechanisms, and an imbalance in sediment organic carbon budgets?, Chem. Rev., 107, 467–485, 2007. 

Carr, M.-E., Friedrichs, M. A., Schmeltz, M., Aita, M. N., Antoine, D., Arrigo, K. R., Asanuma, I., Aumont, O., Barber, R., and Behrenfeld, M.: A comparison of global estimates of marine primary production from ocean color, Deep-Sea Res. Pt. II, 53, 741–770, 2006. 

Cavan, E. L., Trimmer, M., Shelley, F., and Sanders, R.: Remineralization of particulate organic carbon in an ocean oxygen minimum zone, Nat. Commun., 8, 14847,, 2017. 

Chavez, F. P. and Messié, M.: A comparison of eastern boundary upwelling ecosystems, Prog. Oceanogr., 83, 80–96, 2009. 

Cole, S. T., Wortham, C., Kunze, E., and Owens, W. B.: Eddy stirring and horizontal diffusivity from Argo float observations: Geographic and depth variability, Geophys. Res. Lett., 42, 3989–3997, 2015. 

Dale, A. W., Sommer, S., Lomnitz, U., Montes, I., Treude, T., Liebetrau, V., Gier, J., Hensen, C., Dengler, M., Stolpovsky, K., Bryant, L. D., and Wallmann, K.: Organic carbon production, mineralisation and preservation on the Peruvian margin, Biogeosciences, 12, 1537–1559,, 2015. 

DeVries, T. and Weber, T.: The export and fate of organic matter in the ocean: New constraints from combining satellite and oceanographic tracer observations, Global Biogeochem. Cy., 31, 535–555,, 2017. 

Duursma, E. and Boisson, M.: Global oceanic and atmospheric oxygen stability considered in relation to the carbon cycle and to different time scales, Oceanol. Acta, 17, 117–141, 1994. 

Fennel, K., Follows, M., and Falkowski, P. G.: The co-evolution of the nitrogen, carbon and oxygen cycles in the Proterozoic ocean, Am. J. Sci., 305, 526–545, 2005. 

Filippelli, G. M.: The Global Phosphorus Cycle: Past, Present, and Future, Elements, 4, 89–95,, 2008. 

Ganachaud, A. and Wunsch, C.: Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data, Nature, 408, 453,, 2000. 

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate + nitrite, silicate), edited by: Mishonov, A., available at: (last access: November 2019), 2018a. 

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation, edited by: Mishonov, A., available at: (last access: November 2019), 2018b. 

Gruber, N., Friedlingstein, P., Field, C. B., Valentini, R., Heimann, M., Richey, J. E., Lankao, P. R., Schulze, E.-D., and Chen, C.-T. A.: The vulnerability of the carbon cycle in the 21st century: An assessment of carbon-climate-human interactions, Scope-Scientific committee on problems of the environment international council of scientific unions, 62, 45–76, 2004. 

Gruber, N., Frenzel, H., Doney, S. C., Marchesiello, P., McWilliams, J. C., Moisan, J. R., Oram, J. J., Plattner, G.-K., and Stolzenbach, K. D.: Eddy-resolving simulation of plankton ecosystem dynamics in the California Current System, Deep-Sea Res. Pt. I, 53, 1483–1516, 2006. 

Hedges, J. I., Hu, F. S., Devol, A. H., Hartnett, H. E., Tsamakis, E., and Keil, R. G.: Sedimentary organic matter preservation; a test for selective degradation under oxic conditions, Am. J. Sci., 299, 529–555, 1999. 

Henson, S. A., Sanders, R., Madsen, E., Morris, P. J., Le Moigne, F., and Quartly, G. D.: A reduced estimate of the strength of the ocean's biological carbon pump, Geophys. Res. Lett., 38, L04606,, 2011. 

Higgins, M. B., Robinson, R. S., Husson, J. M., Carter, S. J., and Pearson, A.: Dominant eukaryotic export production during ocean anoxic events reflects the importance of recycled NH4+, P. Natl. Acad. Sci USA, 109, 2269–2274,, 2012. 

Inthorn, M., Wagner, T., Scheeder, G., and Zabel, M.: Lateral transport controls distribution, quality, and burial of organic matter along continental slopes in high-productivity areas, Geology, 34, 205–208,, 2006. 

Iversen, M. H. and Ploug, H.: Ballast minerals and the sinking carbon flux in the ocean: carbon-specific respiration rates and sinking velocity of marine snow aggregates, Biogeosciences, 7, 2613–2624,, 2010. 

Iversen, M. H. and Poulsen, L. K.: Coprorhexy, coprophagy, and coprochaly in the copepods Calanus helgolandicus, Pseudocalanus elongatus, and Oithona similis, Mar. Ecol. Prog. Ser., 350, 79–89, 2007. 

Jackson, G. A. and Burd, A. B.: Simulating aggregate dynamics in ocean biogeochemical models, Prog. Oceanogr., 133, 55–65, 2015. 

Karakaş, G., Nowald, N., Schäfer-Neth, C., Iversen, M., Barkmann, W., Fischer, G., Marchesiello, P., and Schlitzer, R.: Impact of particle aggregation on vertical fluxes of organic matter, Prog. Oceanogr., 83, 331–341, 2009. 

Katz, M. E., Finkel, Z. V., Grzebyk, D., Knoll, A. H., and Falkowski, P. G.: Evolutionary Trajectories and Biogeochemical Impacts of Marine Eukaryotic Phytoplankton, Annu. Rev. Ecol. Evol. S., 35, 523–556,, 2004. 

Katz, M. E., Fennel, K., and Falkowski, P. G.: Geochemical and biological consequences of phytoplankton evolution, in: Evolution of primary producers in the sea, Elsevier, 405–430,, 2007. 

Keeling, R. F., Najjar, R. P., Bender, M. L., and Tans, P. P.: What atmospheric oxygen measurements can tell us about the global carbon cycle, Global Biogeochem. Cy., 7, 37–67, 1993. 

Keeling, R. F., Körtzinger, A., and Gruber, N.: Ocean Deoxygenation in a Warming World, Annu. Rev. Mar. Sci., 2, 199–229,, 2010. 

Klausmeier, C. A., Litchman, E., and Levin, S. A.: Phytoplankton growth and stoichiometry under multiple nutrient limitation, Limnol. Oceanogr., 49, 1463–1470, 2004. 

Krumhardt, K. M., Callnan, K., Roache-Johnson, K., Swett, T., Robinson, D., Reistetter, E. N., Saunders, J. K., Rocap, G., and Moore, L. R.: Effects of phosphorus starvation versus limitation on the marine cyanobacterium Prochlorococcus MED4 I: uptake physiology, Environ. Microbiol., 15, 2114–2128,, 2013. 

Kwon, E. Y., Primeau, F., and Sarmiento, J. L.: The impact of remineralization depth on the air–sea carbon balance, Nat. Geosci., 2, 630,, 2009. 

Lam, P. J., Doney, S. C., and Bishop, J. K. B.: The dynamic ocean biological pump: Insights from a global compilation of particulate organic carbon, CaCO3, and opal concentration profiles from the mesopelagic, Global Biogeochem. Cy., 25, GB3009,, 2011. 

Laufkötter, C., John, J. G., Stock, C. A., and Dunne, J. P.: Temperature and oxygen dependence of the remineralization of organic matter, Global Biogeochem. Cy., 31, 1038–1050, 2017. 

Lenton, T. M. and Daines, S. J.: Biogeochemical transformations in the history of the ocean, Annu. Rev. Mar. Sci., 9, 31–58, 2017. 

Lenton, T. M. and Daines, S. J.: The effects of marine eukaryote evolution on phosphorus, carbon and oxygen cycling across the Proterozoic–Phanerozoic transition, Emerging Topics in Life Sciences, 2, 267–278, 2018. 

Lenton, T. M. and Watson, A. J.: Redfield revisited: 1. Regulation of nitrate, phosphate, and oxygen in the ocean, Global Biogeochem. Cy., 14, 225–248, 2000. 

Lenton, T. M., Boyle, R. A., Poulton, S. W., Shields-Zhou, G. A., and Butterfield, N. J.: Co-evolution of eukaryotes and ocean oxygenation in the Neoproterozoic era, Nat. Geosci., 7, 257–265,, 2014. 

Lenton, T. M., Daines, S. J., and Mills, B. J.: COPSE reloaded: an improved model of biogeochemical cycling over Phanerozoic time, Earth-Sci. Rev., 178, 1–28, 2018. 

Lin, S., Litaker, R. W., and Sunda, W. G.: Phosphorus physiological ecology and molecular mechanisms in marine phytoplankton, J. Phycol., 52, 10–36,, 2016. 

Logan, G. A., Hayes, J., Hieshima, G. B., and Summons, R. E.: Terminal Proterozoic reorganization of biogeochemical cycles, Nature, 376, 53–56,, 1995. 

Lomas, M. W., Bonachela, J. A., Levin, S. A., and Martiny, A. C.: Impact of ocean phytoplankton diversity on phosphate uptake, P. Natl. Acad. Sci. USA, 111, 17540–17545,, 2014. 

Lovecchio, E., Gruber, N., Münnich, M., and Lachkar, Z.: On the long-range offshore transport of organic carbon from the Canary Upwelling System to the open North Atlantic, Biogeosciences, 14, 3337–3369,, 2017. 

Lu, W., Ridgwell, A., Thomas, E., Hardisty, D. S., Luo, G., Algeo, T. J., Saltzman, M. R., Gill, B. C., Shen, Y., and Ling, H.-F.: Late inception of a resiliently oxygenated upper ocean, Science, 361, 174–177, 2018. 

Marsay, C. M., Sanders, R. J., Henson, S. A., Pabortsava, K., Achterberg, E. P., and Lampitt, R. S.: Attenuation of sinking particulate organic carbon flux through the mesopelagic ocean, P. Natl. Acad. Sci. USA, 112, 1089–1094, 2015. 

Martiny, A. C., Vrugt, J. A., and Lomas, M. W.: Concentrations and ratios of particulate organic carbon, nitrogen, and phosphorus in the global ocean, Sci. Data, 1, 1–7, 2014. 

McDonnell, A. M. P., and Buesseler, K. O.: Variability in the average sinking velocity of marine particles, Limnol. Oceanogr., 55, 2085–2096,, 2010. 

Meyer, A., Sloyan, B. M., Polzin, K. L., Phillips, H. E., and Bindoff, N. L.: Mixing Variability in the Southern Ocean, J. Phys. Oceanogr., 45, 966–987,, 2015. 

Meyer, K., Ridgwell, A., and Payne, J.: The influence of the biological pump on ocean chemistry: implications for long-term trends in marine redox chemistry, the global carbon cycle, and marine animal ecosystems, Geobiology, 14, 207–219, 2016. 

Moore, J. K., Doney, S. C., and Lindsay, K.: Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model, Global Biogeochem. Cy., 18, GB4028,, 2004. 

Paulmier, A., Kriest, I., and Oschlies, A.: Stoichiometries of remineralisation and denitrification in global biogeochemical ocean models, Biogeosciences, 6, 923–935,, 2009. 

Petit, J.-R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., and Delaygue, G.: Climate and atmospheric history of the past 420 000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436, 1999. 

Ploug, H.: Small-scale oxygen fluxes and remineralization in sinking aggregates, Limnol. Oceanogr., 46, 1624–1631, 2001. 

Polovina, J. J., Mitchum, G. T., and Evans, G. T.: Decadal and basin-scale variation in mixed layer depth and the impact on biological production in the Central and North Pacific, 1960–88, Deep-Sea Res. Pt. I, 42, 1701–1716, 1995. 

Reinhard, C. T., Planavsky, N. J., Gill, B. C., Ozaki, K., Robbins, L. J., Lyons, T. W., Fischer, W. W., Wang, C., Cole, D. B., and Konhauser, K. O.: Evolution of the global phosphorus cycle, Nature, 541, 386–389,, 2017. 

Ridgwell, A.: Evolution of the ocean's “biological pump”, P. Natl. Acad. Sci. USA, 108, 16485–16486, 2011.  

Riley, J., Sanders, R., Marsay, C., Le Moigne, F. A., Achterberg, E. P., and Poulton, A.: The relative contribution of fast and slow sinking particles to ocean carbon export, Global Biogeochem. Cy., 26, 2012. 

Rivkin, R. B. and Legendre, L.: Biogenic Carbon Cycling in the Upper Ocean: Effects of Microbial Respiration, Science, 291, 2398–2400,, 2001. 

Sarmiento, J. L. and Gruber, N.: Ocean Biogeochemical Dynamics, Princeton University Press, Princeton 2006. 

Sharples, J., Middelburg, J. J., Fennel, K., and Jickells, T. D.: What proportion of riverine nutrients reaches the open ocean?, Global Biogeochem. Cy., 31, 39–58,, 2017. 

Slomp, C. P. and Van Cappellen, P.: The global marine phosphorus cycle: sensitivity to oceanic circulation, Biogeosciences, 4, 155–171,, 2007. 

Stukel, M. R. and Ducklow, H. W.: Stirring up the biological pump: Vertical mixing and carbon export in the Southern Ocean, Global Biogeochem. Cy., 31, 1420–1434, 2017. 

Sydeman, W. J., García-Reyes, M., Schoeman, D. S., Rykaczewski, R. R., Thompson, S. A., Black, B. A., and Bograd, S. J.: Climate change and wind intensification in coastal upwelling ecosystems, Science, 345, 77–80,, 2014. 

Tantanasarit, C., Englande, A. J., and Babel, S.: Nitrogen, phosphorus and silicon uptake kinetics by marine diatom Chaetoceros calcitrans under high nutrient concentrations, J. Exp. Mar. Biol. Ecol., 446, 67–75, 2013. 

Van Cappellen, P. and Ingall, E. D.: Redox Stabilization of the Atmosphere and Oceans by Phosphorus-Limited Marine Productivity, Science, 271, 493–496,, 1996. 

Wallmann, K.: Phosphorus imbalance in the global ocean?, Global Biogeochem. Cy., 24, GB4030,, 2010. 

Watson, A. J., Lenton, T. M., and Mills, B. J. W.: Ocean deoxygenation, the global phosphorus cycle and the possibility of human-caused large-scale ocean anoxia, Philos. T. R. Soc. A, 375, 20160318,, 2017. 

Wilson, S. E., Steinberg, D. K., and Buesseler, K. O.: Changes in fecal pellet characteristics with depth as indicators of zooplankton repackaging of particles in the mesopelagic zone of the subtropical and subarctic North Pacific Ocean, Deep-Sea Res. Pt. II, 55, 1636–1647, 2008. 

Wollast, R.: Evaluation and comparison of the global carbon cycle in the coastal zone and in the open ocean, The Sea, 10, 213–252, 1998. 

Yool, A. and Tyrrell, T.: Role of diatoms in regulating the ocean's silicon cycle, Global Biogeochem. Cy., 17, 1103,, 2003. 

Short summary
We present here the newly developed BPOP box model. BPOP is aimed at studying the impact of large-scale changes in the biological pump, i.e. the cycle of production, export and remineralization of the marine organic matter, on the nutrient and oxygen concentrations in the shelf and open ocean. This model has been developed to investigate the global consequences of the evolution of larger and heavier phytoplankton cells but can be applied to a variety of past and future case studies.