The development of deep-ocean anoxia in a comprehensive ocean phosphorus model
TThe development of deep-ocean anoxia in acomprehensive ocean phosphorus model
J. G. Donohue ∗ , B. J. Florio , and A. C. Fowler MACSI, University of Limerick, Limerick, Ireland CSIRO Mineral Resources, Waterford, Western Australia, Australia Complex Systems Group, Department of Mathematics and Statistics,University of Western Australia, Australia OCIAM, University of Oxford, Oxford, UK
Abstract
We analyse a model of the phosphorus cycle in the ocean given by Slompand Van Cappellen (2007, https://doi.org/10.5194/bg-4-155-2007). This modelcontains four distinct oceanic basins and includes relevant parts of the water,carbon and oxygen cycles. We show that the model can essentially be solvedanalytically, and its behaviour completely understood without recourse to nu-merical methods. In particular, we show that, in the model, the carbon andphosphorus concentrations in the different ocean reservoirs are all slaved to theconcentration of soluble reactive phosphorus in the deep ocean, which relaxesto an equilibrium on a time scale of 180,000 y, and we show that the deep oceanis either oxic or anoxic, depending on a critical parameter which we can deter-mine explicitly. Finally, we examine how the value of this critical parameterdepends on the physical parameters contained in the model. The presentedmethodology is based on tools from applied mathematics and can be used toreduce the complexity of other large, biogeochemical models.
Keywords: phosphorus cycle, mathematical model, ocean anoxia event, model re-duction . There are two obvious reasons for wishing to study the phosphorus cycle in the world’soceans. The first is that it is intimately linked to variations in oxygen, carbon and ∗ Corresponding author: [email protected] a r X i v : . [ phy s i c s . a o - ph ] A ug ther elements, both in the atmosphere and in the oceans, and hence also to climate(Van Cappellen and Ingall 1996, Mackenzie et al . 2002). The phosphorus cycle isclosely tied to the biological cycle, particularly in the oceans. While on land eitherphosphorus or nitrogen may be the limiting nutrient, in the ocean it is phosphorusthat is the limiter on geological time scales. This is due to the population of algae(nitrogen fixers) which are able to source nitrogen from the atmosphere (Tyrell 1999).The second reason is that in order to fully understand the effect of anthropogenicalteration of the nitrogen and phosphorus cycle through the use of agricultural fertilis-ers, an understanding of the underlying processes and their time scales of operationis necessary, particularly in view of the impending phosphate crisis (Abelson 1999,Cordell et al . 2009).The phosphorus (or phosphate) cycle has been frequently described (Filipelli 2002,2008, F¨ollmi 1996), but in order to assess and parameterise its effects in the geologicalpast, it is necessary to describe the system using a mathematical model. A numberof such models have been put forward (e. g., Van Capellen and Ingall 1994, Andersonand Sarmiento 1995, Lenton and Watson 2000, Bergman et al . 2004, Tsandev et al .2008, Ozaki et al . 2011), with various applications in mind.One particular application of much recent interest has to do with the occurrenceof ‘oceanic anoxia events’ (OAEs), which have occurred in the geological past, partic-ularly in the Jurassic and Cretaceous periods (Schlanger and Jenkyns 1976, Jenkyns2010). These events are marked in the marine sedimentary record by the occurrenceof organically rich ‘black shales’, and mark periods (of hundreds of thousands of years)during which the deep ocean became anoxic, thus promoting anaerobic digestion andthe production of sulphides and other reduced substances.OAEs are associated with severe changes in climate: warming occurs due to car-bon change in the atmosphere, leading to enhanced precipitation and weathering,hence increased nutrient supply to the oceans, and consequent biomass blooms: thiscauses increased oxygen demand in the upper ocean, and this can lead to deep oceananoxia and eutrophication (Jenkyns 2010). In view of anthropogenic climate change,this raises the question as to whether ocean anoxia is a prospective consequence ofpresent rates of atmospheric carbon increase (Watson 2016). The issue of time scaleprovides another motivation to understand the phosphorus cycle quantitatively, sinceanthropogenic change is drawing down ocean carbonate much more rapidly than itsnormal pre-industrial time scale would suggest (Fowler 2015).It has become increasingly clear that OAEs are frequently associated with theformation of large igneous provinces (LIPs) (Turgeon and Creaser 2008, Sell et al .2014, Percival et al . 2015), and that these may also be associated with increasedweathering (Percival et al . 2016), as well as extinction episodes, which themselvesmight be due to increased upwelling of anoxic water (Jarvis et al . 2008). On theother hand, Van Capellen and Ingall (1994) associate OAEs with a reduced oceancirculation.A number of models have been presented in order to elucidate the circumstances inwhich OAEs occur, for example those of Handoh and Lenton (2003) and Ozaki (2011),and other models have been used for more general purposes (Bergman et al . 2004,Slomp and Van Capellen 2007, Tsandev et al . 2008). A common feature of such models2s their inaccessibility; typically a large number of variables in a number of oceancompartments describe the concentrations of various chemical components, and theseare governed by differential equations which relate changes of the concentrations toreaction terms and inter-compartment fluxes. The complexity of the models is visibleeven in the opacity of their presentation, and their solution is inevitably obtainedthrough numerical simulation. Because of this, it is difficult to interrogate the modelsand virtually impossible to unravel key mechanisms which control the dynamics.The purpose of this paper is to present a methodology, based on tools of appliedmathematics, which can be used to digest such complicated models, and reduce themto a form where their solutions can be obtained cheaply and simply, and the behaviourof the model can be specifically interpreted in terms of the prescribed parameters ofthe model.In particular, we provide an exegesis of the model of Slomp and Van Cappellen(2007), which elaborated the model of Van Cappellen and Ingall (1996) to take ac-count of the difference between continental shelves and the deep ocean. They wereparticularly interested in the effects of ocean mixing on phosphorus burial, and con-sequently on deep ocean anoxia. The numerical results from this model (henceforthcalled the Slomp model) indicate that oxygen concentration and oceanic mixing be-tween these basins significantly affects the phosphorus cycle: in particular, they say:“the simulations show that changes in oceanic circulation may induce marked shiftsin primary productivity and burial of reactive phosphorus between the coastal andopen ocean domains”. Our aim will be to provide explicit parametric interpretationof their results.Our methods, while simple in concept, are sophisticated in practice. They arebased on the ideas of non-dimensionalisation, scaling, and then asymptotic simplifi-cation. As is often the case, the simplifications arise because most of the describingequations act on a faster time scale than the slowest, and thus rate-controlling, equa-tions. This allows us to achieve our goal. In the rest of the paper, the model isdescribed and presented in section 2, and it is then non-dimensionalised in section2.1. The resulting non-dimensional model is incorrectly scaled; we identify the reasonfor this, and correct the problem (by rescaling appropriately). The resulting asymp-totic simplifications are described in 2.2, and lead to the result that all the oceanvariables are slaved to the deep ocean soluble reactive phosphorus, which relaxes toan equilibrium on a time scale of 180,000 y.In section 3, we show that the deep ocean oxygen and reduced substances concen-trations can be determined analytically, and we show that there is a switch from anoxic deep ocean to an anoxic deep ocean at a critical value of one of the dimensionlessparameters. In section 4, we endeavour to unravel the interpretation of our resultsin terms of the physical processes and parameters of the problem; this is the sectionwhere the mathematics-averse should go. Finally we offer our conclusions in section5. We consign much of the algebraic debris to the appendix.3 oastalupwelling DeepSurfaceRiverine EvaporationMixingP. Coastal D. Coastal ExportExport W W W W W F W F W F W F W F W F W F Figure 1: Modelling the oceanic water cycle as four distinct basins with mixing.
The Slomp model divides the ocean into four distinct basins: proximal coastal, dis-tal coastal, surface ocean and deep ocean, having volumes W – W . Volume fluxesbetween the basins are denoted by W F i , i = 1 , , . . . ,
7. The basins and fluxes areshown in figure 1. The fluxes are constructed so that the total inputs and outputsare equal, and the water system is at steady state.The model describes the quantities of phosphorus, carbon and oxygen in the dif-ferent basins. Phosphorus is assumed to be in one of three forms: reactive (SRP),organic particulate (POP), or authigenic calcium phosphate (fish hard parts). Thequantities in each basin are denoted by S i (SRP), P i (POP) and F i (fish P). (Herewe deviate from Slomp and van Cappellen (2007), who allocate P i , i = 1 , , . . . , P F i , i = 1 , , . . . ,
36. The schematic of these processescan be found in figure 2.The carbon cycle is a good deal simpler. It is described by modelling particulateorganic carbon (POC) and is associated with living and detrital biomass. POC maygrow within a water basin depending on phosphorus levels, and additionally there areinter-basin fluxes. The concentration of POC in basin i is denoted by C i . The sourceand flux terms of the carbon system are denoted as CF i , i = 1 , , . . . , W . Thesurface water basins, W , W and W , are assumed to be fully oxic as they are incommunication with the atmosphere. As such, we only model deep ocean oxygenbudget G , which changes in response to to water cycle fluxes, and also aerobic res-piration within W . In an oxygen depleted system, reduced substances like sulphidescan be removed from the system via burial, upwelling or by being oxidised. The con-centration of these reduced substances is denoted by R , and is measured in oxygen4 F P F S P P F P F F P F a P F P F b P F P F P F P F P F P F P F P F b P F P F P S F P F a P F P F P F P F P F P F P F P F P F P F P F P F P F P F S S P F P F Figure 2: A schematic of the model of the oceanic phosphorus cycle. Blue dashedlines show transport via the water cycle, black solid lines indicate biological processeswithin a basin and brown dotted lines represent sedimentary burial.equivalents. There are four source terms OF i for oxygen and five (variously named)for the reduced substances, of which two are flux terms (due to upwelling), and theother three introduce significant nonlinearity. Importantly, the rate of microbial res-piration divides the consumption of deep ocean organic carbon into two components,one of which uses oxygen as the terminal electron acceptor, while the other representsthe use of reduced substances; the split between the two is taken to depend on thedeep ocean oxygen concentration. The full description of the model is given by Slompand van Cappellen (2007), although some of the finer detail is only accessible throughtheir Matlab code.In all, the Slomp model thus consists of eighteen first order differential equationsfor the variables C i , P i , S i , F i , G and R . The quantities X i , X = C, P, F, S, G, R are budgets, i. e., measured in moles, but we prefer to write them as concentrations,thus x i = X i /W i , where W i are the volumes of the basins. Using the definitions ofthe fluxes provided by Slomp and van Cappellen (2007) (and also by Slomp (privatecommunication)), we find that the model takes the form5 c = b s − b c , ˙ c = b c + b s − ( m v + b ) c , ˙ c = b s + ( m v + b ) c − b c , ˙ c = { − h φ ( g ) } ( m v + b ) c + { − h φ ( g ) } b s − b c , ˙ s = b + b p + b f − b s , ˙ p = b s − b p , ˙ f = b s − b f , ˙ s = b s + b p + b f − ( m v + b ) s + vm s , ˙ p = b s + b p − ( m v + b ) p , ˙ f = b s − b f , ˙ s = b p + m vs + ( m v + b ) s − ( m v + b ) s , ˙ p = b s + ( m v + b ) p − ( m v + b ) c − b p , ˙ f = b s − b f , ˙ s = b p + b f + vm s − vm s − b ψ ( g ) , ˙ p = ( m v + b ) (cid:20) R − CP − h χ ( g ) (cid:21) c + b (cid:20) R − CP − h χ ( g ) (cid:21) s − b p , ˙ f = b f − b f , ˙ g = vm ( g s − g ) − b c gK m + g − b Θ( r, g ) , ˙ r = b c (cid:18) − gK m + g (cid:19) − b θ ( r ) − b Θ( r, g ) , (2.1)where we have written g = g and r = r as they have no counterparts in the otherbasins. The coefficients b i and m i are positive constants (with values given in tables4 and 5 of the appendix), R CP is the Redfield ratio of carbon to phosphorus, K m isa Monod constant, g s is the fully oxic surface concentration and v is a dimensionlessmixing parameter. The functions in (2.1) are defined by6 ( r ) = (cid:20) rC RP − (cid:21) + , Θ( r, g ) = 10 − rg,ψ ( g ) = gg for g < g , g ≥ g ,χ ( g ) = .
75 + 0 . gg for g < g , g ≥ g ,φ ( g ) = (cid:18) .
75 + 0 . gg (cid:19) (cid:20) gg + 2371100 (cid:18) − gg (cid:19)(cid:21) − for g < g , g ≥ g , (2.2)where g is a deep oxygen threshold, [ x ] + = max( x, r > g > χ , ψ and φ , only the definitions corresponding to g < g are reported in the article proper.Finally, we note that before the rate law given in Slomp and van Capellen’s equation(4) can be applied, we must first convert r and g from units of moles m − to units ofmoles l − . This leads to the factor of 10 − in (2.2). Our procedure for simplifying the model begins by non-dimensionalising the system.Numerically, the solution approaches a steady state, and our aim is to scale the systemso that the scaled concentration variables are O (1) at the steady state. We first notethat, in Slomp and van Capellen (2007), v = 1 corresponds to the modern, well-mixedocean with values ∼ . v = νv a where v a is a reference value for the anoxic ocean and ν will be our transformed mixingparameter which now ranges between 0 and v a . We set v a = 0 . ν = 1 as a representation of a poorly-mixed ocean with ν = 10 corresponding to themodern, well-mixed ocean.Next, we ensure that the scaled concentrations are O (1) at the steady state byidentifying the largest terms on the right hand sides of system (2.1) and balancingthem. For example, let s = [ s ]¯ s , where [ s ] denotes the scale, and ¯ s is the newdimensionless variable. For some equations, the scaling argument is straightforward;7onsider for example (2.1) , which becomes[ c ][ t ] ˙¯ c = b [ s ]¯ s − b [ c ]¯ c , (2.3)where [ t ] is the chosen time scale. A balance of the two terms on the right hand sidegives 0 = b [ s ] − b [ c ] , (2.4)which relates the scales of [ s ] and [ c ]. For equations with more than two terms, theresults of a numerical simulation are used to infer the largest two terms. One mustbe careful in some situations where a cyclic definition of scales is found. For example,consider (2.1) and (2.1) , where taking the largest two terms gives0 = b [ p ] − b [ s ] , b [ s ] − b [ p ] , (2.5)for which the only solution is [ p ] = [ s ] = 0. To resolve this conundrum, we consideralso the next largest term. In this particular case, it is the constant riverine input,giving 0 = b + b [ p ] − b [ s ] , b [ s ] − b [ p ] , (2.6)which has a non-trivial solution. Physically, this occurs because a large amount ofphosphorus is cycled between SRP and POP phases compared to the net input andoutput. Two further instances of this cyclicity occur in the choice of scales for s and s . It is perhaps easier to see how scales are chosen by restricting ourselves to themost obvious balances. These are r ∼ C RP , − rg = Θ ∼ m v a g s b ,s ∼ b c b ∼ b p b ∼ b f b ,s ∼ b c b ∼ b f b ∼ b p b ,s ∼ b c b ∼ b c b ∼ b p b ∼ b f b ∼ b R CP p b ,s ∼ b p m v a , f ∼ b f b , t ∼ m v a ∼ ,
000 y , (2.7)with the time scale being chosen as the longest time scale of any of the equations,which leads to the consequence that all of the time derivatives (bar that of the slowestequation) will be multiplied by parameters less than one (and in fact much less thanone). The values associated with these scales are given in (A.1) of the appendix. Theresulting scaled system is given by 8 ˙ c = s − c ,ε ˙ c = δ c + s − ε νc − c ,ε ˙ c = s + ε νc + ε c − c ,ε ˙ c = ( ε + ε ν ) c − ( ε + ε ν ) c φ ( g ) + s − ε s φ ( g ) − c ,ε ˙ s = λ + λ p + ε f − s ,ε ˙ p = s − p ,ε ˙ f = s − f ,ε ˙ s = ε s + λ p + ε f + ε νs − ε νs − s ,ε ˙ p = s + δ p − ε νp − p ,ε ˙ f = s − f ,ε ˙ s = λ p + δ νs + δ s + λ νs − ε νs − s ,ε ˙ p = s + ε p + ε νp − ε c − ε νc − p ,ε ˙ f = s − f , ˙ s = p + ε f + ε νs − s ν − λ ψ ( g ) ,ε ˙ p = ( ε + ε ν ) c − ( ε + ε ν ) c χ ( g ) + s − ε s χ ( g ) − p ,ε ˙ f = f − f ,ε ˙ g = ν (1 − ε g ) − ε c gλ + g − rg,ε ˙ r = ε c (cid:18) − gλ + g (cid:19) − [ r − + − δ rg, (2.8)where we have omitted the overbars for convenience. The functions in these equationsare defined by φ ( g ) = (cid:18) .
75 + 0 . gg ∗ (cid:19) (cid:26) gg ∗ + 2371100 (cid:18) − gg ∗ (cid:19)(cid:27) − for g < g ∗ , g ≥ g ∗ ,ψ ( g ) = gg ∗ for g < g ∗ , g ≥ g ∗ ,χ ( g ) = .
75 + 0 . gg ∗ for g < g ∗ , g ≥ g ∗ , (2.9)where g ∗ = g / [ g ]. 9he dimensionless coefficients are defined in (A.3)–(A.5) in the appendix. Theyare divided into three sets; parameters denoted λ i are of O (1); parameters denoted δ i are small ∼ .
1, but not very small; and parameters ε i are ‘very small’, in practice < .
1. There is some fuzziness at the crossover, for example the parameters ε , ε , ε , ε , ε , ε , ε , ε , ε and ε could all have been taken as δ s.The scales in (2.7) give fifteen of the eighteen scales necessary, and it can be seenthat of the equations, no precise balance has been applied in the equations for s , s and s . As explained above, the scale for s is chosen by solving (2.6); this isequivalent to choosing λ = 1 − λ . (2.10)In a similar manner, the scales for s and s are chosen by defining λ = 1 − ε ,λ = 1 − δ . (2.11)This then completes the choice of scaling of the model. To determine if the scalingis appropriate for a poorly-mixed ocean, we now compute the dimensionless steadystate solution with ν = 1; denoting these values with an asterisk, these are found tobe g ∗ = 0 . , c ∗ = 1 . ,c ∗ = 26 . , c ∗ = 783 . ,c ∗ = 774 . , s ∗ = 1 . ,p ∗ = 1 . , f ∗ = 1 . ,s ∗ = 26 . , p ∗ = 26 . ,f ∗ = 26 . , s ∗ = 779 . ,p ∗ = 783 . , f ∗ = 779 . ,s ∗ = 906 . , p ∗ = 782 . ,f ∗ = 779 . , r ∗ = 1 . . (2.12)We might expect that the steady state values would be O (1), but clearly this is not thecase. It seems that there is a magnifying factor of about 30 from basin 1 to basin 2,and then 30 from basin 2 to basin 3. There is some subtle effect here, which needs tobe elucidated. There are two key scales: [ s ] and [ s ]; every other scale can be relatedback to these. We focus on the steady state solutions of the differential equationsfor s and s . Substituting in the other variables and neglecting small terms, we candeduce ( ε ν + ε − ε ) s − ε s = ( λ δ + ε ) s − λ ε ψ ( g ) , − (( λ ε + λ ) ν + δ ) s + ( ε ν + δ − δ ) s = λ δ ε s ν − δ λ ψ ( g ) . (2.13)The coefficients of s and s on the left hand side of these equations form a 2 × ≈ . ν = 1. This explains why the10ystem is sensitive to inaccuracies. When ν = 1, the values of the diagonals of thismatrix are ε + ε − ε ≈ .
029 and ε + δ − δ ≈ . s ] = 1 ε + ε − ε , [¯ s ] = 1( ε + ε − ε )( ε + δ − δ ) . (2.14)Thus, from the original dimensionless variables, we now define ¯ s = [¯ s ]ˆ s , ¯ s = [¯ s ]ˆ s ,and from these we can deduce the rescaling of all the other variables other than r , g and those in basin 1, which are unaltered, just as in (2.7). The rescaled system isnow found to be (in terms of the hatted variables, but again we drop the hats) ε ˙ c = s − c ,ε ˙ c = ε c + s − ε νc − c ,ε ˙ c = s + ε νc + ε c − c ,ε ˙ c = ( ε + ε ν ) c − ( ε + ε ν ) c φ ( g ) + s − ε s φ ( g ) − c ,ε ˙ s = λ + λ p + ε f − s ,ε ˙ p = s − p ,ε ˙ f = s − f ,ε ˙ s = ε s + λ p + ε f + ε νs − ε νs − s ,ε ˙ p = s + ε p − ε νp − p ,ε ˙ f = s − f ,ε ˙ s = λ p + δ νs + ε s + ε νs − ε νs − s ,ε ˙ p = s + ε p + ε νp − ε c − ε νc − p ,ε ˙ f = s − f , ˙ s = p + ε f + ε νs − s ν − ε ψ ( g ) ,ε ˙ p = ( ε + ε ν ) c − ( ε + ε ν ) c χ ( g ) + s − ε s χ ( g ) − p ,ε ˙ f = f − f ,ε ˙ g = ν (1 − ε g ) − λ c gλ + g − rg,ε ˙ r = δ c (cid:18) − gλ + g (cid:19) − [ r − + − δ rg. (2.15)The new dimensionless coefficients are defined in (A.6). With ν = 1, the steady-state11olution in these rescaled variables is given by g ∗ = 0 . , c ∗ = 1 . ,c ∗ = 0 . , c ∗ = 0 . ,c ∗ = 0 . , s ∗ = 1 . ,p ∗ = 1 . , f ∗ = 1 . ,s ∗ = 0 . , p ∗ = 0 . ,f ∗ = 0 . , s ∗ = 0 . ,p ∗ = 0 . , f ∗ = 0 . ,s ∗ = 0 . , p ∗ = 0 . ,f ∗ = 0 . , r ∗ = 1 . . (2.16)As they are now all O (1), it shows that the current scaling is adequate for a poorly-mixed ocean. Inspecting (2.15), it is clear that on a rapid time scale ( ∼ ε i ), s , c , p , f → c ≈ s ≈ p ≈ f , but the degeneracy betweenthe s and p equations leaves their value indeterminate. As is usual in this situation,the missing information is obtained by eliminating the large term; we add the s and p equations, and this leads to (bearing in mind the basin 1 and basin 2 equalities)( ε + ε ) ˙ s = ε + ε + ε νs − ( ε − ε + ν ( ε + ε )) s , (2.17)suggesting a slower evolution of the basin 2 variables. Similarly, the basin 3 concen-trations all rapidly equilibrate, but there is degeneracy in the s and p equations,and adding these yields( ε + ε ) ˙ s = δ νs + ( ε + ε − ε + ( ε + ε − ε ) ν ) s − ( δ + ε ν ) s . (2.18)Finally, the basin 4 variables f , c , p → s rapidly, and thus the slow s equation is˙ s ≈ ( λ + ε ν ) s − νs , (2.19)where λ = 1 + ε ∼ . . (2.20)Thus, we can write the s and s equations in the form ε ˙ s = δ + νs − ( λ + λ ν ) s ,ε ˙ s = νs + ( δ + δ ν ) s − ( λ + δ ν ) s , (2.21)12here ε = ε + ε ε ∼ . , ε = ε + ε δ ∼ . × − ,δ = ε + ε ε ∼ . , δ = ε + ε − ε δ ∼ . ,δ = ε + ε − ε δ ∼ . , δ = ε δ ∼ . ,λ = ε − ε ε ∼ . , λ = ε + ε ε ∼ . ,λ = δ δ ∼ . . (2.22)We have broken our rule about the size of δ s and ε s, but it is necessary to retain theapparently small terms in δ , δ and δ . Evidently the s and s equations are stillrelatively fast, and clearly both of them relax to an equilibrium approximately givenby s ≈ δ + s νλ + λ ν , s ≈ (( δ + λ ) ν + δ + λ ) νs + δ ( δ + δ ν )( λ + δ ν )( λ + λ ν ) , (2.23)following which s relaxes to its equilibrium s ≈ ( ε ν + λ ) δ ( δ ν + δ )( δ λ − ( δ + λ ) ε ) ν + ε ν + ( λ ( λ − λ ) − λ δ ) ν ∼ .
943 (2.24)with ε = ( λ − λ ) λ − ε ( δ + λ ) − δ λ + δ λ ∼ . × − . (2.25)This relaxation occurs on a time scale t ∼ ( δ ν + λ )( λ + λ ν )( δ λ − ( δ + λ ) ε ) ν + ε ν + ( λ ( λ − λ ) − λ δ ) ν ∼ . , (2.26)corresponding to 180,000 y, much longer than our original time scale. Numericalverification of these analytical estimates is given in figure 3 where the four SRPvariables are used as exemplars. Although the equations for p , s and c are coupled to g in (2.15), the coupling is weakand can be ignored, so that the transport part of the model operates independentlyfrom oxygen and reduced substances in the deep ocean. The model equations for r and g are given by the last pair in (2.15), and depend on the transport only through c ≈ s , which is given by (2.23), and varies on a slow time scale. Thus ε ˙ g = ν (1 − ε g ) − λ s gλ + g − rg,ε ˙ r = δ s (cid:18) − gλ + g (cid:19) − [ r − + − δ rg. (3.1)13 − − − − − − − − − − − log time (years)log time (years) log time (years)log time (years) s ( mM ) s ( mM ) s ( mM ) s ( mM ) Figure 3: Equilibration of soluble reactive phosphorus concentrations in each oceanicbasin. All variables are presented in dimensional form and time has been logarith-mically transformed in order to clearly illustrate the various time scales of interest.We have set the mixing parameter ν = 1 to represent a poorly-mixed ocean. Theinitial value of all eighteen system variables was set to be 1 and (2.15) was solvednumerically. 14 g r Figure 4: The approximate slow manifold g = G ( r ), or g nullcline, (3.2), usingparameter values from the Slomp model and ν set to 1 to represent a poorly-mixedocean. In the lower curve, s = 0 .
782 whereas in the upper curve, s = 0 . ε ∼ − whereas ε ∼ − and therefore the g equation relaxes first to anequilibrium in which r + ε ν = νg − λ s λ + g . (3.2)This defines g as a decreasing function G ( r ), with G (0) being finite or very large (as ε ∼ − ) depending on whether λ s > ν or < ν respectively; figure 4 shows twotypical examples, one with s = 0 .
782 (corresponding to the steady state in (2.16))and one using a much smaller value of s .Following the relaxation of g to its quasi-equilibrium G ( r ), r evolves (still relativelyrapidly) via the ε ˙ r equation, which can be written, using (3.2), in the approximateform ε ˙ r = δ ( λ s − ν + ε νg ) λ + (cid:18) δ λ − δ (cid:19) rg − [ r − + . (3.3)Figure 5 plots ε ˙ r as a function of r for the two values of s used in figure 4. We seethat for the normal value s = 0 . r andthus g are O (1), and because of our choice of scales the deep ocean is anoxic, but thatfor s < νλ ≈ . r collapses to zero, and the oxygen level increases dramatically.This sharp transition is due to the apparent parametric accident that δ λ − δ = 0,according to the values in (A.4) and (A.6). It seems unlikely such a coincidence wouldoccur, but in fact, working our way through the definitions of the parameters in the15 ǫ ˙ r r Figure 5: ε ˙ r as a function of r given by (3.4) using parameter values from theSlomp model and ν set to 1 to represent a poorly-mixed ocean. In the upper curve, s = 0 .
782 whereas in the lower curve, s = 0 . δ λ δ = 1 . Ultimately, this is due to the equal coefficients b in the rates of aerobic and anaerobicrespiration in (2.1). The model thus takes the very simple form in the anoxic case λ s > ν : ε ˙ r = δ (cid:18) s − νλ + ε νgλ (cid:19) − [ r − + . (3.4)Anoxic equilibrium is obtained in a time scale t ∼ ε ∼ − , corresponding to about30 y. What if λ s < ν ? It is then necessary to rescale the variables as r ∼ ε , g ∼ ε , (3.5)and (3.1) now takes the approximate form (since ε (cid:28) ε ε ˙ g = ν − λ s − νg − rg,ε ε ˙ r = − δ rg ; (3.6)thus r → t ∼ ε /ε ∼ g → − λ s ν , and in dimensional terms, 0 . − λ s ν ) mM.16 .2 Numerical verification Section 3.1 provided a description of the dynamical behaviour of the oxygen subsys-tem in oxic and anoxic conditions as well as characterising the transition between theoxic and anoxic states. We will assess the accuracy of these predictions through nu-merical solutions of (2.15). Slomp and Van Cappellen (2007) showed that the mixingparameter could be varied to induce switches between oxic and anoxic conditions.Thus, in figure 6, we plot steady-state values of g and r as a function of the mixingparameter ν . Note that here, and in the remainder of this section, we have reverted todimensional variables for ease of interpretation. Examining the numerical solutions,we note that at a critical value of ν ≈ . r falls abruptly from 0 .
03 mM to near zerowhile g begins to increase rapidly. Thus, the sudden shift in (equilibrium) redox statepredicted once a critical parameter value has been exceeded (see section 3) appearsto be borne out by the numerical solutions. r ( m M ) ν g ( m M ) ν Figure 6: Steady-state values of r and g as a function of the dimensionless oceanicmixing parameter ν . Solutions of (2.15) were obtained numerically before variableswere converted to their dimensional forms.To verify that we have successfully captured the mechanism behind this abruptchange, we use our numerical output to plot ˙ r as a function of r . We carry out thisexercise both sides of the apparent discontinuity with the results shown in figure 7.A small change in ν brings about a drastic shift in the position of the ˙ r curve andhence a large change in the equilibrium value of r . It is instructive to compare thesecurves with the dimensionless equivalents in figure 5 which have assumed ν = 1. Itappears that the mixing parameter is sufficiently high in figure 7 that the ε νg termin (3.4) is no longer negligible. This has the effect of converting the flat piece of ˙ r to amonotonically decreasing function of r . Nonetheless, the relationship between r and˙ r at low r values is relatively insensitive, facilitating the large shift in steady-stateconcentration as the mixing parameter moves through a critical threshold.Finally, using ν = 1 and ν = 4 . −3 ˙ r ( m M / y ) r (mM) Figure 7: ˙ r as a function of r plotted at ν = 4 . ν = 4 . r and all other system variables are at their numericallyobtained steady-state values.30 y. Meanwhile, in section 3.1, we predicted that a well-mixed deep ocean wouldrecover its oxygen levels in a time scale of approximately 3000 y. In order to assess thevalidity of these estimate, we set all other variables to steady-state and plot numericalsolutions for g and r with ν = 1 (see figure 8(i)) and ν = 4 . The analysis in section 2.2 suggests that the chemical components in the differentcompartments rapidly (here meaning (cid:28) ,
000 y) come to an approximate equilibrium,where the values are determined in terms of the deep ocean reactive phosphorus s ,which however evolves over a much longer time scale ∼ ,
000 y to an eventualequilibrium given by (2.24). The surface ocean reactive phosphorus s follows thesame slow evolution, being determined by (2.23). During this slow evolution of s ,the deep ocean will rapidly (30 y) become anoxic if λ s > ν , whereas if λ s < ν itbecomes oxic, slightly less rapidly (3,000 y).Put simply, our analysis suggest that when s is too large, the deep ocean willbecome anoxic. Assuming the ocean is generally near steady-state conditions, wehave a statement involving the equilibrium value of s . We compute this equilibriumvalue both numerically and using our analytical approximation and plot λ s − ν asa function of ν in figure 9. By comparison with the Slomp article’s “degree of anox-icity” measure, we observe that this quantity successfully captures the deep ocean’s18
75 1501.41.61.822.22.4 x 10 − − r ( mM ) g ( mM ) r ( mM ) g ( mM ) time (years) time (years) (ii)(i) Figure 8: Equilibration of oxygen ( g ) and reduced substances ( r ) concentrations attwo different mixing rates ( ν ). In both cases, a steady state for the overall system ofdifferential equations is first found numerically with the associated concentrations ofoxygen and reduced substances then perturbed to 80% of their steady-state values.In (i), ν = 1 and the ocean is poorly mixed whereas in (ii), ν = 4 . ν ≈
4. This model therefore has thecapacity to explain ocean anoxia events, depending on the assumed parameters ofthe problem. It is thus important to unravel what all these complicated parametercombinations mean in terms of the dimensional parameters of the physical system.While λ is independent of the mixing parameter ν , the value of s depends on ν aswell as other system parameters. This functional dependence is not known exactly.However, by using the approximate form of the s steady-state value, we can express λ s − ν as an explicit function of the model’s dimensional parameter set.Unfortunately, although the analysis is simple, the dependence of the criticalparameter on the physical inputs is non-trivial in the extreme, to the extent that inthe appendix we give an algorithm to compute λ s approx − ν , and in the electronicsupplementary material provide a code which does this (see Online Resource 1). Thefully expanded expression for λ s approx − ν depends on 48 dimensional parameters.Here, we focus on the influence of b , the riverine input of SRP. Slomp and VanCappellen’s (2007) numerical exploration revealed that anoxia would occur if thepresent ocean’s circulation rate was reduced by 50% ( ν = 5 in our notation) whilethe supply of reactive phosphorus from the continents was simultaneously boosted by20%. They suggested that such an increase could be caused by coastal erosion linkedto sealevel rise.Setting ν = 5 and leaving all other parameters at their previously assigned values,we plot λ s approx − ν as a function of b . Figure 10 demonstrates that λ s approx − ν switches from negative to positive as b is increased with the crossover occurringwhen b is 6% higher than its baseline value. Numerical study (not shown) revealsthat the threshold actually lies at a value of b that is 12% higher than the baselinevalue (i.e., between our prediction and the value used by Slomp and Van Cappellen).19 ano x i c t h r e s ho l d s λ s approx − νλ s num − ν DOA − ν Figure 9: Anoxia onset metrics as a function of the mixing parameter ν . The quantity λ s − ν is presented in two forms, one using a numerical estimate of the s steady-state value ( λ s num − ν ) and one using an analytical approximation ( λ s approx − ν ).Negative values of the metrics correspond to oxic conditions. The degree of anoxicity(DOA) quantity, defined as 1 − gg ∗ in the original Slomp paper, is shown as a reference.Thus, the quantity λ s approx − ν appears to be able to predict changes in ocean oxygenstatus, whether they be linked to circulation or other factors. In this article, we have systematically analysed the model of the phosphorus cyclein the ocean given by Slomp and Van Cappellen (2007). Through careful scaling ofthe Slomp model, we identified a large number of negligible steady-state fluxes. Wealso isolated distinct time scales associated with system equilibration. By exploitingthese two factors, we were able to effectively decouple the subsystem of oxygen andreduced substances from the carbon-phosphorus cycle.While soluble reactive phosphorus acts as an (effectively static) input to the oxy-gen subsystem, the contribution of oxygen to the cycling of carbon and phosphoruscan be safely ignored. In particular, this means that a range of nonlinear, non-smoothfunctions used to model redox dependence in the burial of sorbed P, particulate or-ganic P and particulate organic carbon can be excised without affecting our qualita-tive findings. From a starting of point of eighteen nonlinear equations, we separatelyanalysed a set of sixteen (approximately) linear equations which govern carbon andphosphorus transport and a pair of equations which explain the chemistry of the oxicdeep ocean, the chemistry of the anoxic deep ocean and the nature of the transitionbetween the two.Having partitioned the system into two parts, we can elucidate the nature of the20 .5 2.6 2.7 2.8x 10 −3 −0.25−0.1500.150.25 λ s a pp r o x − ν b (riverine input) Figure 10: Anoxia onset metric as a function of the riverine input parameter b wherewe have set ν = 5. Positive values of λ s − ν correspond to anoxic conditions andnegative values to oxic.transition between oxic and anoxic oceans. A small change in system parametersproduces abrupt, almost discontinuous, switches in the equilibrium concentrations ofoxygen and reduced substances. We link this sensitivity in the model to the functionalform prescribed for microbial respiration. Allison and Martiny (2008) refer to this kindof microbial model as a “black box” with “microorganisms buried within equationstructure as kinetic constants and response functions”. Our analysis highlights theneed to compare the predictions of such studies with those of models that explicitlyincorporate microbial biomass, in order to enhance understanding of how anoxiaoccurs.With the nature of transition to anoxia determined, we sought to determine thesystem parameters responsible for driving such a transition. Unfortunately, due tothe scope of the original model (containing 69 parameters), the critical controllingparameter defies succinct characterisation. However, by focusing on a small subset ofthe system parameters (i.e., mixing rate and riverine input), we demonstrated thatone can accurately predict the outcome of changes in the the rate of a given process.While our focus here was on providing this kind of proof-of-concept, future workcould entail analytic study of the expression (A.7) in the appendix. In particular,one can explicitly determine whether the ocean’s oxygen status is affected by vari-ation (or covariation) in a few parameters of interest. More generally, we suggestthat this article demonstrates the viability of adopting a systematic, mathematicalapproach in studying the behaviour of large, biogeochemical models. The deductionof parameterised steady-state concentrations and equilibration time scales, as we havepresented here, is generally beyond the reach of a purely computational approach tobiogeochemistry. 21 cknowledgements We acknowledge the support of the Mathematics Applications Consortium for Sci-ence and Industry ( ) funded by the Science Foundation Irelandmathematics grant 12/1A/1683. This publication has emanated from research con-ducted with the financial support of Science Foundation Ireland under Grant Number:SFI/09/IN.1/I2645. No new data was used in this study.
Appendix
Non-dimensionalisation
We non-dimensionalise the model as described before (2.8). In the model as presentedby Slomp and Van Capellen (2007) the chemical quantities are measured in Tmol andthe reservoir volumes in Tm . We have written the model equations (2.1) in termsof concentrations, which thus have the units mol m − = mM (millimolar). The set ofscales is given by (in units of mM)[ r ] = C RP ≈ × − , [ s ] = b d ≈ . × − , [ c ] = b b [ s ] ≈ . × − , [ c ] = b b [ s ] ≈ . × − , [ s ] = d [ s ] ≈ . × − , [ s ] = d [ s ] ≈ . × − , [ c ] = b b [ s ] ≈ . × − , [ c ] = b b [ s ] ≈ . × − , [ p ] = b b [ s ] ≈ . × − , [ f ] = b b [ s ] ≈ . × − , [ p ] = b b [ s ] ≈ . × − , [ f ] = b b [ s ] ≈ . × − , [ p ] = b b [ s ] ≈ . × − , [ f ] = b b [ s ] ≈ . × − , [ s ] = b m v a [ p ] ≈ . × − , [ p ] = b R CP b [ s ] ≈ . × − , [ f ] = b b [ f ] ≈ . × − , [ g ] = 10 m v a g s b [ r ] ≈ . × − . (A.1)where d = b − b b b , d = b b b b − b b , d = b b b b − b b . (A.2)22he dimensionless parameters of the resulting model (2.8) are given by the followingthree sets. First the O (1) parameters λ i : λ = b b [ s ] ≈ . , λ = b [ p ] b [ s ] ≈ . , λ = b [ p ] b [ s ] ≈ . ,λ = b [ p ] b [ s ] ≈ . , λ = b [ s ] m v a ≈ . , λ = K m [ g ] ≈ . ,λ = m v a [ s ] b [ s ] ≈ . . (A.3)Next, the small, but not very small, parameters δ i : δ = m v a [ s ] b [ s ] ≈ . , δ = b [ s ] b [ s ] ≈ . , δ = b [ r ][ g ]10 b ≈ . ,δ = b [ c ] b [ c ] ≈ . δ = b [ p ] b [ p ] ≈ . . (A.4)23inally, the very small parameters ε i : ε = b [ f ][ s ] b ≈ . × − , ε = b [ s ] b [ s ] ≈ . × − ,ε = b [ f ] b [ s ] ≈ . × − , ε = m v a [ s ] b [ s ] ≈ . × − ,ε = b [ f ][ s ] m v a ≈ . × − , ε = b [ p ][ p ] b ≈ . × − ,ε = b [ c ][ p ] b ≈ . × − , ε = h b [ c ]237[ p ] b ≈ . × − ,ε = h b [ s ]237[ p ] b ≈ . × − , ε = b [ c ][ c ] b ≈ . × − ,ε = b [ c ][ c ] b ≈ . × − , ε = b [ c ] h [ c ] b ≈ . × − ,ε = b h [ s ][ c ] b ≈ . × − , ε = [ g ] g s ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = m [ g ] v a m g s ≈ . × − , ε = b [ c ][ p ] b R CP ≈ . × − ,ε = m v a b ≈ . × − ε = m v a b ≈ . × − ε = m v a b ≈ . × − ε = m v a [ r ] b ≈ . × − ε = b [ c ] b ≈ . × − , ε = b [ c ] m v a g s ≈ . × − ,ε = m v a [ c ][ p ] b R CP ≈ . × − , ε = m v a b ≈ . × − ,ε = m v a b ≈ . × − , ε = m v a b ≈ . × − ,ε = h m [ c ] v a p ] b ≈ . × − , ε = m [ c ] h v a [ c ] b ≈ . × − , (A.5)24 = m v a b ≈ . × − , ε = [ s ][ s ] ≈ . × − ,ε = m v a [ p ][ p ] b ≈ . × − , ε = m v a [ c ][ p ] b ≈ . × − ,ε = m v a [ c ][ c ] b ≈ . × − , ε = m v a [ c ][ c ] b ≈ . × − . When the rescaling introduced before (2.15) is done, the new dimensionless pa-rameters are given by λ = ε [¯ s ] ≈ . , δ = ε [¯ s ] ≈ . ,ε = δ [¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = δ [¯ s ] ≈ . × − , ε = δ [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = λ [¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = ε [¯ s ][¯ s ] ≈ . × − , ε = ε [¯ s ][¯ s ] ≈ . × − ,ε = λ [¯ s ][¯ s ] ≈ . × − . (A.6) Dimensional parameters
In this section, we present the definitions and values of all dimensional constantsassociated with the model (2.1). Before we do so, we first document the parametersof the original Slomp model. We begin with parameters governing the water cycle.As shown in table 1, the fluxes corresponding to river input (
W F ), ocean upwelling25abel Flux Definition W F River input
W k W F Proximal to distal =
W F W F Distal to surface =
W F + W F W F Ocean downwelling =
W F + W F W F Ocean upwelling v o W k W F Coastal upwelling v c W k W F Evaporation =
W F Table 1: Definition of water fluxes in the Slomp model. The values of the constants
W k i , v o and v c are given in table 2.( W F ) and coastal upwelling ( W F ) are defined empirically, via constants that wewill refer to as W k , W k and W k respectively. Changes in circulation are modelledby multiplying the oceanic and coastal upwelling constants by the non-dimensionalparameters v o and v c , respectively. The remaining four fluxes in the oceanic circula-tion system then arise by imposing conservation of stationary water volume withineach basin.The values assigned to the water cycle parameters, along with all other parametersof the Slomp model are listed in tables 2 and 3. These parameters combine to producethe dimensional constants used in (A.1) and the model (2.1). The definitions of thesedimensional constants are given in tables 4 and 5.26arameter Value Units W . × m W . × m W . × m W . × m W k . × m /yW k . × m /yW k . × m /yv o vv c vCk . × y − Ck . y − Ck . × − Ck Ck . y − Ck . y − Ck . / (560 .
25 + 4 . Ck Ck . × − y − Ck . × − y − Ck (496 . / (3600 + 28 . Ck . × − y − R CP . × R CO / k redox . × mM y − k prec × − mM y − [RS] P k × mol/y P k . × − y − P k . × − P k P k a × − P k b × − y − P k c P k .
43 y − P k P k P k . × − y − Table 2: Definition of the parameters of the Slomp model (continued below)27arameter Value Units
P k . × − P k P k a × − P k b × − y − P k c P k .
18 y − P k . / (0 .
044 + 5 . P k P k . × − y − P k × − P k P k × − y − P k . × − y − P k × − y − P k . × mol/ yP k . × − P k . / . P k m . × − y − ( W k + W k ) /W b . × − y − Ck /R CO b . × mM y − k redox b . × y − (1 − Ck ) × Ck × R CP b .
97 y − ( Ck W + Ck W k ) /W b . × − y − (cid:0) (1 − Ck ) Ck W k (cid:1) /W b . × y − (1 − Ck ) Ck R CP b .
21 y − ( Ck + Ck W k ) /W m . × − y − W k Ck /W b . × y − (1 − Ck ) Ck R CP b . × − y − − Ck Ck W k /W + Ck W k /W m . × − y − (1 − Ck ) W k Ck /W b . × − y − Ck b . × − y − Ck Ck W k /W m . × − y − Ck W k Ck /W b . × − y − ( Ck Ck R CP W ) /W b . × − y − Ck b .
01 y − (1 − P k ) P k b × − y − P k b b . × y − (cid:0) ( Ck + P k ) W + P k W k (cid:1) /W b . × y − Ck (cid:0) − P k a − ( R CP P k Ck / (cid:1) b .
46 y − ( P k W + P k W k ) /W b . × − y − P k a Ck b × − y − P k b b . × − y − ( P k W k ) /W b .
18 y − (1 − P k ) P k b × − y − P k b b .
07 y − Ck + P k + P k W k /W m . × − y − P k W k /W b .
05 y − Ck (1 − P k − P k a ) b . × − y − P k W k (1 − P k ) /W b .
19 y − P k + P k W k /W m . × − y − P k W k /W b . × − y − P k a Ck b × − y − P k b b . × − y − P k b . × − y − P k W k /W m . × − y − P k W k /W b . × − y − Ck m . × − y − ( W k + W k ) /W Table 4: Definition of constants (continued below).29onstant Value Units Definition b . × − y − Ck (1 − P k − Ck ) b . × − y − P k W k /W m . × − y − P k W k /W b . × − y − ( Ck Ck W k ) / ( R CP W )) m . × − y − ( Ck Ck W k ) / ( R CP W )) b . × − y − P k b . × − y − P k Ck b × − y − P k b . × − y − P k (1 − P k ) b × − y − P k b . × − y − Ck Ck W k /W m . × − y − Ck W k Ck /W b . × − y − ( Ck Ck W R CP ) /W b . × − y − P k b . × − y − (cid:0) P k (1 − P k ) W (cid:1) /W b × − y − P k b × − mM y − k prec b . × − mM y − P k /W m . × − y − W k /W m . × − y − W k /W m . × − y − ( W k + W k ) /W b . × − mM y − P k /W h . × − Ck C RP g .
17 mM [O ] t =0 g s .
325 mM k O2–surf R CP K m × − mM v a . W i is thevolume of basin i and W k i , Ck i and P k i are rate constants in the water, carbon andphosphorus cycles respectively. 30 noxia parameter computation Although our analysis of the model provided us with an extremely simple result (theoxygen status of the deep ocean depending on the sign of λ s − ν ), the determinationof the critical parameters involved in the transition to anoxia is convoluted in the ex-treme. Therefore here we provide a path to compute them, and in the supplementarymaterial we provide a Matlab code to compute them directly (see Online Resource1), given the original input parameters of the Slomp model (those listed in tables 2and 3).From these, tables 4 and 5 provide definitions of all b i and m i , h , R CP , C RP , K m , g and g s . From these, (A.1) provides sequential definitions for all the scales [ r ], [ s ],etc., where additionally (2.10) and (2.11) have been used; d , d and d are definedin (A.2). From these, (A.3) defines λ ,. . . λ , λ and λ , (A.4) defines δ ,. . . δ , (A.5)defines ε ,. . . ε , ε , ε ,. . . ε , ε ,. . . ε , ε ,. . . ε , ε , ε . . . ε and ε . . . ε .We then use (2.14) to define ¯ s and ¯ s , after which (A.6) defines λ , δ , ε ,. . . ε and ε ,. . . ε . Finally, we recover λ from (2.20), ε from (2.25) and (2.22) gives ε , ε , δ , . . . δ and λ , . . . λ .The net result of these transformations is that, in dimensional terms, λ s approx − ν can be expressed as b b X b + X b + b b b b b b b b m b b R CP b ( ν ( m b b − m b b ) v a − b b b + b b b ) b − v a m ( X b + ( m b b b b b ( m b νv a − b b ) R CP + X b ) b b ) g s b b ( b b − b b ) , (A.7) X = − b ( b (( − m m + m ( m − m )) b + b m ( m − m )) b + m m b b b ) b b ν b ( b b − b b ) g s m R CP m b b v a − b ν R ( b b − b b ) g s m v a − ν R v a − b (( b b b b b + b b b b b ) b − b b b b b b ) b b b b b m R CP b b ,X = ν m m m b m b b b b b R CP g s ( b b − b b ) v a + b b ν ( b b − b b ) g s ((( b ( m − m ) b + b b m ) b m b + b b b b m m ) b R CP + m m b b b b b ) m v a − ν (( − ( b b − b b ) g s m (( m b b b b b + b b b ( m b b − m b b )) b R CP + b b b b ( m b b − m b b )) b + m m b b b b b b b b b R CP ) b + m m b b b b b b b R CP g s ( b b − b b )) v a − m b b R CP b b b b b b b b b ,X = − b ((( m b νv a b ( m νv a + b ) b + ( ν ( m − m ) v a − b ) b b b b ) b − b b b (( m νv a + b ) b − b b )) m b b + b b ( m b b ( m νv a + b ) b + m b b b b ) b b ) R CP − b ( m b b ( m νv a + b ) b + m b b b b ) b b b b ,X = b ((( ν (( m m + m ( m − m )) b + b m ( m − m )) v a + ν (( b m − m b + b ( m − m )) b − b b m − b b ( m − m )) v a − b ( b b − b b )) b + νm b b v a ( m νv a + b )) m b b + b (( ν (( m m + m m ) b + m b m ) v a + ( m b + m b ) b − m b b ) b + m b b ( m νv a + b )) b b ) R CP + (( ν (( m m + m m ) b + m b m ) v a +( m b + m b ) b − m b b ) b + m b b ( m νv a + b )) b b b b , R = (( b (( b m (( m b − b m + b ( m − m )) b + b b m ) b − b (( m m + m m ) b + m b m ) b b ) b R CP − b b b b (( m m + m m ) b + m b m )) b + b b b m ( b ( m b b b + m b b b ) R CP + m b b b b )) b − m b b R CP b b b b b b ( m − m )) b − m b b b b b b b R CP ( m b + m b ) ,R = ( − b (((( − m b b b b + b b b ( m b + m b )) b R CP + b b b b ( m b + m b )) b b − b b b m ( R CP b b b b + b b b b )) b − b ( b ( m b b b − m b b b ) R CP + m b b b b ) b b b )( b b − b b ) g s m b + b b b b b m R CP (( m b b b b + m b b b b ) b − m b b b b b ) b b ) b − m m b b R CP b b b b b b g s ( b b − b b )( b b − b b ) . Thus, our efforts to write an explicit formula for λ s approx − ν lead to extremely com-plicated formulae having no apparent simplification; it thus appears that the simplecontrolling parameters of the solution depend in a very complicated way on manyof the physically prescribed parameters, and this dependence needs to be elucidatedcomputationally (see Online Resource 1). References
Abelson, P. H. 1999 A potential phosphate crisis. Science (5,410), 2015.Allison, S. D. and J. B. Martiny 2008 Resistance, resilience, and redundancy in mi-crobial communities. PNAS (1), 11512–11519.Anderson, L. A. and J. L. Sarmiento 1995 Global ocean phosphate and oxygen sim-ulations. Global Biogeochem. Cycles (4), 621–636.Bergman, N. M., T. M. Lenton and A. J. Watson 2004 COPSE: a new model ofbiogeochemical cycling over Phanerozoic time. Amer. J. Sci. , 397–437.Cordell, D., J.-O. Drangert and S. White 2009 The story of phosphorus: global foodsecurity and food for thought. Global Environmental Change , 292–305.Filipelli, G. ,M. 2002 The global phosphorus cycle. In: Kohn, M., J. Rakovanand J. Hughes (eds.) Phosphates: geochemical, geobiological, and materialsimportance. Revs. Mineral. Geochem. , pp. 391–425.Filipelli, G. ,M. 2008 The global phosphorus cycle: past, present, and future. Ele-ments , 89–95. 33¨ollmi, K. B. 1996 The phosphorus cycle, phosphogenesis and marine phosphate-richdeposits. Earth-Sci. Revs. , 55–124.Fowler A. C. 2015 A simple thousand year prognosis for ocean and atmosphere car-bon change. Pure Appl. Geophys. (1), 49–56.Handoh, I. C. and T. M. Lenton 2003 Periodic mid-Cretaceous oceanic anoxic eventslinked by oscillations of the phosphorus and oxygen biogeochemical cycles.Global Biogeochem. Cycles (4), 1092.Jarvis, I., G. A. Carson, M. K. E. Cooper, M. B. Hart, P. N. Leary, B. A. Tocher,D. Horne and A. Rosenfeld 1988 Microfossil assemblages and the Cenomanian-Turonian (late Cretaceous) oceanic anoxic event. Cretaceous Res. , 3–103.Jenkyns, H. C. 2010 Geochemistry of oceanic anoxic events. Geochemistry Geo-physics Geosystems (3), Q03004.Lenton, T. M. and A. J. Watson 2000 Redfield revisited: 2. What regulates theoxygen content of the atmosphere? Global Biogeochem. Cycles , 249–268.Mackenzie, F. T., L. M. Vera and A. Lerman 2002 Century-scale nitrogen and phos-phorus controls of the carbon cycle. Chem. Geol. , 13–32.Ozaki, K., S. Tajima and E. Tajika 2011 Conditions required for oceanic anoxia/euxinia:constraints from a one-dimensional ocean biogeochemical cycle model. EarthPlanet. Sci. Letts. , 270–279.Percival, L. M. E., A. S. Cohen, M. K. Davies, A. J. Dickson, S. P. Hesselbo, H. C.Jenkyns, M. J. Leng, T. A. Mather, M. S. Storm and W. Xu 2016 Osmium iso-tope evidence for two pulses of increased continental weathering linked to EarlyJurassic volcanism and climate change. Geology (9), 759–762.Percival, L. M. E., M. L. I. Witt, T. A. Mather, M. Hermoso, H. C. Jenkyns, S. P.Hesselbo, A. H. Al-Suwaidi, M. S. Storm, W. Xu and M. Ruhl 2015 Globally en-hanced mercury deposition during the end-Pliensbachian extinction and Toar-cian OAE: a link to the Karoo?Ferrar Large Igneous Province. Earth Planet.Sci. Letts. , 267–280.Schlanger, S. O. and H. C. Jenkyns 1976 Cretaceous oceanic anoxic events: causesand consequences. Geologie en Mijnbouw (3-4), 179–184.Sell, B., M. Ovtcharova, J. Guex, A. Bartolini, F. Jourdan, J. E. Spangenberg, J.-C. Vicente and U. Schaltegger 2014 Evaluating the temporal link between theKaroo LIP and climatic-biologic events of the Toarcian Stage with high-precisionU-Pb geochronology. Earth Planet. Sci. Letts. , 48–56.Slomp, C. P. and P. van Capellen 2007 The global marine phosphorus cycle: sensi-tivity to oceanic circulation. Biogeosci. , 155–171.34sandev, I., C. P. Slomp and P. van Capellen 2008 Glacial-interglacial variations inmarine phosphorus cycling: implications for ocean productivity. Global Bio-geochem. Cycles , GB4004.Turgeon, S. C. and R. A. Creaser 2008 Cretaceous oceanic anoxic event 2 triggeredby a massive magmatic episode. Nature , 323–326.Tyrrell, T. 1999 The relative influence of nitrogen and phosphorus on oceanic pri-mary production. Nature , 525–531.Van Cappellen, P. and E. D. Ingall 1994 Benthic phosphorus regeneration, net pri-mary production, and ocean anoxia: a model of the coupled marine biogeo-chemical cycles of carbon and phosphorus. Paleoceanography (5), 677–692.Van Cappellen, P. and E. D. Ingall 1996 Redox stabilization of the atmosphere andoceans by phosphorus-limited marine productivity. Science (5,248), 493–496.Watson, A. J. 2016 Oceans on the edge of anoxia. Science354