The development of deep-ocean anoxia in a comprehensive ocean phosphorus model

We analyse a model of the phosphorus cycle in the ocean given by Slomp and Van Cappellen (Biogeosciences 4:155–171, 2007. 10.5194/bg-4-155-2007). This model contains four distinct oceanic boxes and includes relevant parts of the water, carbon and oxygen cycles. We show that the model can essentially be solved analytically, and its behaviour completely understood without recourse to numerical methods. In particular, we show that, in the model, the carbon and phosphorus concentrations in the different ocean reservoirs are all slaved to the concentration of soluble reactive phosphorus in the deep ocean, which relaxes to an equilibrium on a time scale of 180,000 y, and we show that the deep ocean is either oxic or anoxic, depending on a critical parameter which we can determine explicitly. Finally, we examine how the value of this critical parameter depends on the physical parameters contained in the model. The presented methodology is based on tools from applied mathematics and can be used to reduce the complexity of other large, biogeochemical models. Supplementary Information The online version contains supplementary material available at 10.1007/s13137-023-00221-0.


Introduction
There are two obvious reasons for wishing to study the phosphorus cycle in the world's oceans.The first is that it is intimately linked to variations in oxygen, carbon and other elements, both in the atmosphere and in the oceans, and hence also to climate (Van Cappellen andIngall 1996, Mackenzie et al. 2002).The phosphorus cycle is closely tied to the biological cycle, particularly in the oceans.While on land either phosphorus or nitrogen may be the limiting nutrient, in the ocean it is phosphorus that is believed to be 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 anthropogenic alteration of the nitrogen and phosphorus cycle through the use of agricultural fertilisers, an understanding of the underlying processes and their time scales of operation is 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öllmi 1996), but in order to assess and parameterise its effects in the geological past, it is necessary to describe the system using a mathematical model.A number of such models have been put forward (e. g., Van Capellen and Ingall 1994, Anderson and 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 occurrence of 'oceanic anoxia events' (OAEs), which have occurred in the geological past, particularly in the Jurassic and Cretaceous periods (Schlanger andJenkyns 1976, Jenkyns 2010).These events are marked in the marine sedimentary record by the occurrence of organically rich 'black shales', and mark periods (of hundreds of thousands of years) during which the deep ocean became anoxic, thus promoting anaerobic digestion and the production of sulphides and other reduced substances.
It has become increasingly clear that OAEs are frequently associated with the formation 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 increased weathering (Percival et al. 2016), as well as extinction episodes, which themselves might be due to increased upwelling of anoxic water (Jarvis et al. 2008).
OAEs are also associated with severe changes in climate: warming occurs due to carbon change in the atmosphere, leading to enhanced precipitation and weathering, hence increased nutrient supply to the oceans, and consequent biomass blooms: this causes increased oxygen demand in the upper ocean, and this can lead to deep ocean anoxia (Jenkyns 2010).Eutrophic conditions in the surface ocean may be further enhanced by redox-dependent release of phosphorus from anoxic sediments (Van Capellen and Ingall 1994).In view of anthropogenic climate change, this raises the question as to whether ocean anoxia is a prospective consequence of present rates of atmospheric carbon increase (Watson 2016).On the other hand, Niemeyer et al. (2017) suggest that the positive benthic P-release feedback may be mitigated by the configuration of the modern ocean, preventing a full-scale OAE.
It is clear that the mechanisms through which OAEs are sustained are controversial (Beil et al. 2020) with evidence often generated through the simulation of detailed numerical box models, for example those of Handoh and Lenton (2003), Slomp and Van Capellen (2007) and Wallmann et al. (2019).Thus, there is a need to enhance understanding of how these models produce a prediction, rather than allowing them to become black boxes (Maeda et al. 2021).Unfortunately, a common feature of such models is their inaccessibility; typically a large number of variables in a number of oceanic 'boxes' describe the concentrations of various chemical components, and these are governed by differential equations which relate changes of the concentrations to reaction terms and inter-box fluxes.The complexity of the models is visible even in the opacity of their presentation, and their solution is inevitably obtained through numerical simulation.Because of this, it is difficult to interrogate the models and virtually impossible to unravel key mechanisms which control the dynamics.
The purpose of this paper is to present a methodology, based on tools of applied mathematics, which can be used to digest such complicated models, and reduce them to a form where their solutions can be obtained cheaply and simply, and the behaviour of the model can be specifically interpreted in terms of the prescribed parameters of the 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 account of the difference between continental shelves and the deep ocean.They were particularly interested in the effects of ocean mixing on phosphorus burial, and consequently on deep ocean anoxia.The numerical results from this model (henceforth called the Slomp model) indicate that oxygen concentration and mixing between boxes significantly affects the phosphorus cycle: in particular, they say: "the simulations show that changes in oceanic circulation may induce marked shifts in primary productivity and burial of reactive phosphorus between the coastal and open ocean domains".Our aim will be to provide explicit parametric interpretation of their results.
Our methods, while simple in concept, are sophisticated in practice.They are based on the ideas of non-dimensionalisation, scaling, and then asymptotic simplification.As is often the case, the simplifications arise because most of the describing equations act on a faster time scale than the slowest, and thus rate-controlling, equations.This allows us to achieve our goal.In the rest of the paper, the model is described and presented in section 2, and it is then non-dimensionalised in section 2.1.The resulting non-dimensional model is incorrectly scaled; we identify the reason for this, and correct the problem (by rescaling appropriately).The resulting asymptotic simplifications are described in 3.1, and lead to the result that all the ocean variables are slaved to the deep ocean soluble reactive phosphorus, which relaxes to an equilibrium on a time scale of 180,000 y.
In section 3.2, we show that the deep ocean oxygen and reduced substances concentrations can be determined analytically, and we show that there is a switch from an oxic deep ocean to an anoxic deep ocean at a critical value of one of the dimensionless parameters.In section 4, we endeavour to unravel the interpretation of our results in terms of the physical processes and parameters of the problem; this is the section where the mathematics-averse should go.Finally we offer our conclusions in section 5. We consign much of the algebraic debris to the appendix.The model describes the quantities of phosphorus, carbon and oxygen in the different oceanic boxes.Phosphorus is assumed to be in one of three forms: reactive (SRP), organic particulate (POP), or authigenic calcium phosphate (fish hard parts).The quantities in each box are denoted by S i (SRP), P i (POP) and F i (fish P). (Here we deviate from Slomp and van Cappellen (2007), who allocate P i , i = 1, 2, . . ., 12 to these variables.)The phosphorus budgets are altered either by reactive processes within an oceanic box, or by travelling from one box to another.

Label Flux Definition
The carbon cycle is a good deal simpler.It is described by modelling particulate organic carbon (POC) and is associated with living and detrital biomass.POC may grow within an oceanic box depending on phosphorus levels, and additionally there are inter-box fluxes.The concentration of POC in box i is denoted by C i .
The modelling of the oxygen system is assumed to be important only for the deep ocean, W 4 .The surface-level boxes, W 1 , W 2 and W 3 , are assumed to be fully oxic as they are in communication with the atmosphere.As such, we only model deep ocean oxygen budget G 4 , which changes in response to water-cycle fluxes, and also aerobic respiration within W 4 .In an oxygen depleted system, reduced substances like sulphides can be removed from the system via burial, upwelling or by being oxidised.The concentration of these reduced substances is denoted by R 4 , and is measured in oxygen equivalents.Importantly, the rate of microbial respiration divides the consumption of deep ocean organic carbon into two components, one of which uses oxygen as the terminal electron acceptor, while the other represents the use of reduced substances; the split between the two is taken to depend on the deep ocean oxygen concentration.The full description of the model is given by Slomp and van Cappellen (2007), although some of the finer detail is only accessible through their Matlab code.
In all, the Slomp model thus consists of eighteen first order differential equations for the variables C i , P i , S i , F i , G 4 and R 4 .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 boxes.The fluxes of the water, phosphorus, carbon and oxygen cycles are associated with a set of parameters denoted W k i , P k i , Ck i and Ok i respectively.Their description and values are given in tables 2 and 3 of the appendix.The conversion of moles to concentrations produces a transformed version of this parameter set which is described in tables 4 and 5 of the appendix.We find that the converted model takes the form where we have written g 4 = g and r 4 = r as they have no counterparts in the other boxes.The coefficients a i and b i are positive constants related to the phosphorus and carbon cycle respectively, whereas the coefficients m i are positive constants that do not fit neatly into either of the former categories.Their values are given in tables 4 and 5 of the appendix.Furthermore, R CP is the Redfield ratio of carbon to phosphorus with (C/P ) oxic and (C/P ) anoxic the ratios of carbon to phosphorus for sedimentary organic matter buried under oxic and anoxic conditions, respectively.Finally, k redox controls the reduced-substances reoxidation rate, k prec controls the removal of reduced substances via precipitation, K m is a Monod constant, g s is the fully oxic surface concentration and v is a dimensionless mixing parameter.The functions in (2.1) are defined by for g < g 0 , where g 0 is a deep oxygen threshold, [x] + = max(x, 0), and we have assumed that r > 0 and g > 0 in the definition of Θ.These functions correspond to the flux terms given in equations ( 3)-( 7) of Slomp and van Capellen (2007).However, for χ, ψ and ϕ, only the definitions corresponding to g < g 0 are reported in the article proper.Note that the factor of 0.25 in the definitions of χ and ϕ arises from the assumption that anoxia may reduce the burial flux of p 4 by up to 25%.Finally, we note that before the rate law (describing reoxidation of reduced substances) given in Slomp and van Capellen's equation ( 4) can be applied, we must first convert r and g from units of moles m −3 to units of moles l −1 .This leads to the factor of 10 −6 in (2.2).

Non-dimensionalisation of the model
Our procedure for simplifying the model begins by non-dimensionalising the system.Numerically, using the parameter values estimated by Slomp and van Capellen (2007), the solution approaches a steady state.Our aim is to scale the system so that the scaled concentration variables are O(1) at this steady state.We first note that in Slomp and van Capellen's model, the mixing parameter v was taken to be = 1 for a well-mixed ocean, but was lower for poorly mixed anoxic oceans, with values v ∼ v a = 0.1.We use this anoxic value in our choice of scales below, but because later in section 3.2 we also consider the case of a well-mixed ocean, it is useful to retain the dimensionless parameter in order to facilitate the possibility of adjusting the mixing parameter in a straightforward manner.
Next, we ensure that the scaled concentrations are O(1) at the steady state by identifying the largest terms on the right hand sides of system (2.1) and balancing them.For example, let s 1 = [s 1 ]s 1 , where [s 1 ] denotes the scale, and s1 is the new dimensionless variable.For some equations, the scaling argument is straightforward; consider for example (2.1) 1 , which becomes where [t] is the chosen time scale.A balance of the two terms on the right hand side gives which relates the scales of [s 1 ] and [c 1 ].For equations with more than two terms, the results of a numerical simulation are used to infer the largest two terms.One must be careful in some situations where a cyclic definition of scales is found.For example, consider (2.1) 5 and (2.1) 6 , where taking the largest two terms gives for which the only solution is [p 1 ] = [s 1 ] = 0. To resolve this conundrum, we consider also the next largest terms.In this particular case, it is the constant riverine input, giving which has a non-trivial solution.Physically, this occurs because a large amount of phosphorus is cycled between SRP and POP phases compared to the net input and output.Two further instances of this cyclicity occur in the choice of scales for s 2 and s 3 .It is perhaps easier to see how scales are chosen by restricting ourselves to the most obvious balances.These are 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 slowest equation) will be multiplied by parameters less than one (and in fact much less than one).The values associated with these scales are given in (A.1) of the appendix.The resulting scaled system is given by where we have omitted the overbars for convenience.The functions in these equations are defined by 1 for g ≥ g * 0 , (2.10) The dimensionless coefficients are defined in (A.3)-(A.5) in the appendix.They are divided into three sets; parameters denoted λ i are of O(1); parameters denoted δ i are small ∼ 0.1, but not very small; and parameters ε i are 'very small', in practice < 0.1.There is some fuzziness at the crossover, for example the parameters ε 6 , ε 8 , ε 13 , ε 14 , ε 33 , ε 99 , ε 107 , ε 110 , ε 112 and ε 113 could all have been taken as δs.
The scales in (2.8) give fifteen of the eighteen scales necessary, and it can be seen that of the equations, no precise balance has been applied in the equations for s 1 , s 2 and s 3 .As explained above, the scale for s 1 is chosen by solving (2.7); this is equivalent to choosing λ 2 = 1 − λ 1 . (2.11) In a similar manner, the scales for s 2 and s 3 are chosen by defining (2.12) This then completes the choice of scaling of the model.To determine if the scaling is appropriate for a poorly-mixed ocean, we now compute the dimensionless steady state solution with ν = 1; denoting these values with an asterisk, these are found to be We might expect that the steady state values would be O(1), but clearly this is not the case.Inspecting the sixteen carbon and phosphorus variables, it seems that there is a magnifying factor of about 30 from box 1 to the corresponding box 2 variable, and then 30 from box 2 to the corresponding box 3 variable.There is some subtle effect here, which needs to be elucidated.There are two key scales: [s 2 ] and [s 3 ]; every other scale can be related back to these.We focus on the steady state solutions of the differential equations for s 2 and s 3 .Substituting in the other variables and neglecting small terms, we can deduce The coefficients of s 2 and s 3 on the left hand side of these equations form a 2 × 2 matrix which has a small determinant (≈ 0.0007) when ν = 1.This explains why the system is sensitive to inaccuracies.When ν = 1, the values of the diagonals of this matrix are ε 101 + ε 2 − ε 3 ≈ 0.029 and ε 106 + δ 2 − δ 1 ≈ 0.034.In order to accommodate the fact that these numbers are very small, it is appropriate to rescale the variables.We do this by defining rescaling parameters . (2.15) Thus, from the original dimensionless variables, we now define s2 = [s 2 ]ŝ 2 , s3 = [s 3 ]ŝ 3 , and from these we can deduce the rescaling of all the other variables other than r, g and those in box 1, which are unaltered, just as in (2.8).The rescaled system is now found to be (in terms of the hatted variables, but again we drop the hats) The new dimensionless coefficients are defined in (A.6).With ν = 1, the steady-state solution in these rescaled variables is given by (2.17) As they are now all O(1), it shows that the current scaling is adequate for a poorlymixed ocean.

Model reduction
In this section, we study the dynamics of the scaled Slomp model.It is important to note that we will assume that the Slomp model has been well parameterised.Specifically, we will not allow for the possibility that their estimates of the system parameters differ from the 'true' values by an order of magnitude or more.On this basis, a number of simplifications to the model can be made, as we will see in sections 3.1 and 3.2.It is also important to note that, in the analysis that follows, we will assume that the (dimensionless) initial conditions are O(1) or equivalently that all variables are within an order of magnitude of their equilibrium values, as given by (2.17).Through numerical investigation, it is apparent that this is the only stable steady-state solution associated with these parameter values.However, it is, of course, possible that there are additional steady-state solutions at other parameter values, as similar box models have been shown to exhibit bistable equilibria (Goldblatt et al. 2006) and sustained oscillations (Handoh andLenton 2003, Wallmann et al. 2019).If such dynamics were to occur in the Slomp model, they could be traced back to one or more of the five nonlinear equations in (2.16).Thus, our approach does not rule out the detection of more complex dynamics, though a search is not a focal part of the analysis.

Simplification of the carbon-phosphorus model
Inspecting (2.16), it is clear that on a rapid time scale ( Similarly, in box 2, we rapidly have c 2 ≈ s 2 ≈ p 2 ≈ f 2 , but the degeneracy between the s 2 and p 2 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 2 and p 2 equations, and this leads to (bearing in mind the box 1 and box 2 equalities) suggesting a slower evolution of the box 2 variables.Similarly, the box 3 concentrations all rapidly equilibrate, but there is degeneracy in the s 3 and p 3 equations, and adding these yields Finally, the box 4 variables f 4 , c 4 , p 4 → s 3 rapidly, and thus the slow where Thus, we can write the s 2 and s 3 equations in the form where We have broken our rule about the size of δs and εs, but it is necessary to retain the apparently small terms in δ 7 , δ 8 and δ 9 .Evidently the s 2 and s 3 equations are still relatively fast, and clearly both of them relax to an equilibrium approximately given by following which s 4 relaxes to its equilibrium  This relaxation occurs on a time scale corresponding to 180,000 y, much longer than our original time scale.Numerical verification of these analytical estimates is given in figure 2 where the four SRP variables are used as exemplars.

Oxygen dynamics
Although the equations for p 4 , s 4 and c 4 are coupled to g in (2.16), the coupling is weak and can be ignored, so that the carbon-phosphorus part of the model operates independently from oxygen and reduced substances in the deep ocean.The model equations for r and g are given by the last pair in (2.16), and depend on the carbon and phosphorus equations only through c 4 ≈ s 3 , which is given by (3.7), and varies on a slow time scale.Thus Now ε 32 ∼ 10 −4 whereas ε 37 ∼ 10 −2 and therefore the g equation relaxes first to an equilibrium in which This defines g as a decreasing function G(r), with G(0) being finite or very large (as ε 19 ∼ 10 −4 ) depending on whether λ 6 s 3 > ν or < ν respectively; figure 3 shows two typical examples, one with s 3 = 0.782 (corresponding to the steady state in (2.17)) and one using a much smaller value of s 3 .Following the relaxation of g to its quasi-equilibrium G(r), r evolves (still relatively rapidly) via the ε 37 ṙ equation, which can be written, using (3.12), in the approximate form Figure 4 plots ε 37 ṙ as a function of r for the two values of s 3 used in figure 3. We see that for the normal value s 3 = 0.782, there is a stable steady state in which r and thus g are O(1), and because of our choice of scales the deep ocean is anoxic.However, for s 3 < ν λ 6 ≈ 0.33, r collapses to zero, and the oxygen level increases dramatically.
This sharp transition is due to the apparent parametric accident that δ 25 λ 6 − δ 3 = 0, according to the values in (A.4) and (A.6).It seems unlikely such a coincidence would occur, but in fact, working our way through the definitions of the parameters in the appendix, we do find that δ 3 λ 6 δ 25 = 1. (3.14) Ultimately, this is due to the equal coefficients b 1 in the rates of aerobic and anaerobic respiration in (2.1).The model thus takes the very simple form in the anoxic case λ 6 s 3 > ν: Anoxic equilibrium is obtained in a time scale t ∼ ε 37 ∼ 10 −2 , corresponding to about 30 y.

Numerical verification
We have provided a description of the dynamical behaviour of the oxygen subsystem in oxic and anoxic conditions as well as characterising the transition between the oxic and anoxic states.We will now assess the accuracy of these predictions through numerical solutions of (2.16).Slomp and Van Cappellen (2007) showed that the mixing parameter could be varied to induce switches between oxic and anoxic conditions.Thus, in figure 5, we plot steady-state values of g and r as a function of the mixing parameter ν.Note that here, and in the remainder of this section, we have reverted to dimensional variables for ease of interpretation.Examining the numerical solutions, we note that at a critical value of ν ≈ 4.14, r falls abruptly from 0.03 mM to near zero while g begins to increase rapidly.Thus, the sudden shift in (equilibrium) redox state predicted once a critical parameter value has been exceeded (see section 3.2) appears to be borne out by the numerical solutions.In view of the preceding discussion, one might expect to see a jump in g at the same critical value of ν ≈ 4.14 in figure 5.This is masked by the fact that the jump in g corresponding to the jump in r is on the anoxic oxygen scale ∼ 3 × 10 −5 mM, and thus not visible in figure 5. Further, we can see from (3.15) that the jump in r occurs when λ 6 s 3 ≈ ν, and that from (3.12), the anoxic oxygen is g ≈ λ 11 ν λ 6 s 3 − ν when r ≈ 0. So when r jumps down, the anoxic-scaled g increases rapidly, and this is indicated by the rapid rise in g (on the oxic scale) in figure 5, which is proportional to ν − λ 6 s 3 , as can be seen from (3.17).
To verify that we have successfully captured the mechanism behind this abrupt change, we use our numerical output to plot ṙ as a function of r.We carry out this exercise on both sides of the apparent discontinuity with the results shown in figure 6.A small change in ν brings about a drastic shift in the position of the ṙ curve and hence a large change in the equilibrium value of r.It is instructive to compare these curves with the dimensionless equivalents in figure 4 which have assumed ν = 1.It appears that the mixing parameter is sufficiently high in figure 6 that the ε 19 νg term in (3.13) is no longer negligible.This has the effect of converting the flat piece of ṙ to a monotonically decreasing function of r.Nonetheless, the relationship between r and ṙ at low r values is relatively insensitive, facilitating the large shift in steady-state concentration as the mixing parameter moves through a critical threshold.Finally, using ν = 1 and ν = 4.5 to represent anoxic and oxic oceans respectively, we examine the dynamics of the oxygen sub-system.We recall that, in section 3.2, we predicted that anoxic equilibrium would be obtained in a time scale of approximately 30 y.Meanwhile, in section 3.2.1,we predicted that a well-mixed deep ocean would recover its oxygen levels in a time scale of approximately 3000 y.In order to assess the validity of these estimates, we set all other variables to steady-state and plot numerical solutions for g and r with ν = 1 (see figure 7(i)) and ν = 4.5 (see figure 7(ii)).In both cases, we observe strong agreement between the numerical output and (ii) (i) Figure 7: Equilibration of oxygen (g) and reduced substances (r) concentrations at two different mixing rates (ν).In both cases, a steady state for the overall system of differential equations is first found numerically with the associated concentrations of oxygen 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.5 and the ocean more closely resembles the present-day configuration.In both cases, we solve (2.16) numerically and then convert all variables to their dimensional forms.our analytical predictions.

Discussion of transition to anoxia
The analysis in section 3.1 suggests that the chemical components in the different oceanic boxes rapidly (here meaning ≪ 3,000 y) come to an approximate equilibrium, where the values are determined in terms of the deep ocean reactive phosphorus s 4 , which however evolves over a much longer time scale ∼ 180,000 y to an eventual equilibrium given by (3.8).The surface ocean reactive phosphorus s 3 follows the same slow evolution, being determined by (3.7).During this slow evolution of s 3 , the deep ocean will rapidly (30 y) become anoxic if λ 6 s 3 > ν, whereas if λ 6 s 3 < ν it becomes oxic, slightly less rapidly (3,000 y).
In section 3.2, we analysed the mechanisms responsible for shifts between anoxic and oxic deep oceans in the model.Starting with a poorly-mixed ocean, we observed that the processes of reoxidation of reduced substances and aerobic respiration appear in both the differential equation for oxygen and the differential equation for reduced substances.Equilibrium of the oxygen equation implies that the losses of oxygen to these two processes are effectively cancelled out by the net supply of oxygen from the surface ocean.This, in turn, means that the remaining two terms in the reducedsubstances equation must balance at equilibrium.One of these terms is a constant (or weakly decreasing) input of reduced substances.The other term corresponds to the removal of reduced substances via precipitation (assumed to be followed by burial in sediment).However, this precipitation, modelled by k prec θ(r) in the dimensional system, is activated only when the concentration of reduced substances exceeds a The quantity λ 6 s 3 −ν is presented in two forms, one using a numerical estimate of the s 3 steady-state value (λ 6 s num 3 − ν) and one using an analytical approximation (λ 6 s approx 3 − ν).Negative values of the metrics correspond to oxic conditions.The degree of anoxicity (DOA) quantity, defined as 1 − g g * 0 in the original Slomp paper, is shown as a reference.Note that oxygen levels, though extremely small at low values of the mixing parameter, remain nonzero and therefore, the DOA never equals one.prescribed threshold value.This non-smooth feature of the model produces a kink in the relationship between ṙ and r (see figure 4 and figure 6).The presence of this kink means that large changes in equilibrium concentrations can be brought about by small changes in system parameters.
In simple biogeochemical terms, our analysis suggests that when surface ocean reactive phosphorus s 3 is too large, the deep ocean will become anoxic.Assuming the ocean is generally near steady-state conditions, we have a statement involving the equilibrium value of s 3 .We compute this equilibrium value both numerically and using our analytical approximation (given by (3.7)) and plot λ 6 s 3 − ν as a function of ν in figure 8.We note that the upper-limit value of ν = 10 corresponds to the modern, well-mixed ocean.By comparison with the Slomp article's 'degree of anoxicity' measure, we observe that this quantity successfully captures the deep ocean's transition from an oxic to an anoxic state at ν ≈ 4.This model therefore has the capacity to explain ocean anoxia events, depending on the assumed parameters of the problem.It is thus important to unravel what all these complicated parameter combinations mean in terms of the dimensional parameters of the physical system.While λ 6 is independent of the mixing parameter ν, the value of s 3 depends on ν as well as other system parameters.This functional dependence is not known exactly.However, by using the approximate form of the s 3 steady-state value, we can express λ 6 s 3 − ν as an explicit function of the model's dimensional parameter set.
Unfortunately, although the analysis is simple, the dependence of the critical parameter on the physical inputs is non-trivial in the extreme, to the extent that in the appendix we give an algorithm to compute λ 6 s approx 3 − ν, and in the electronic supplementary material provide a code which does this (see Online Resource 1).The fully expanded expression for λ 6 s approx 3 − ν depends on 48 dimensional parameters.Here, we focus on the influence of a 53 , the riverine input of SRP.Slomp and Van Cappellen's (2007) numerical exploration revealed that anoxia would occur if the present ocean's circulation rate was reduced by 50% (ν = 5 in our notation) while the supply of reactive phosphorus from the continents was simultaneously boosted by 20%.They suggested that such an increase could be caused by coastal erosion linked to sea level rise.
Setting ν = 5 and leaving all other parameters at their previously assigned values, we plot λ 6 s approx 3 − ν as a function of a 53 .Figure 9 demonstrates that λ 6 s approx 3 − ν switches from negative to positive as a 53 is increased with the crossover occurring when a 53 is 6% higher than its baseline value.Numerical study (not shown) reveals that the threshold actually lies at a value of a 53 that is 12% higher than the baseline value (i.e., between our prediction and the value used by Slomp and Van Cappellen).Thus, the quantity λ 6 s approx 3 − ν appears to be able to predict changes in ocean oxygen status, whether they be linked to circulation or other factors.
It is important to emphasise that we focused on changes in ν and a 53 only because Slomp and Van Cappellen's (2007) numerical exploration examined these two factors.The purpose of producing an algebraic expression for λ 6 s approx 3 − ν was to understand how the transition to anoxia depends on system parameters, more broadly.This 'All-At-A-Time' approach to sensitivity analysis (Pianosi et al. 2016) allows us to assess the robustness of the model output to the modelling assumptions while avoiding the need to produce a high volume of model runs and then visually compare model predictions.

Conclusions
In this article, we have systematically analysed the model of the phosphorus cycle in the ocean given by Slomp and Van Cappellen (2007).Through careful scaling of the Slomp model, we identified a large number of negligible steady-state fluxes.We also isolated distinct time scales associated with system equilibration.By exploiting these two factors, we were able to effectively decouple the subsystem of oxygen and reduced substances from the carbon-phosphorus cycle.
While soluble reactive phosphorus acts as an (effectively static) input to the oxygen subsystem, the contribution of oxygen to the cycling of carbon and phosphorus can be safely ignored.In particular, this means that a range of nonlinear, non-smooth functions used to model redox dependence in the burial of sorbed P, particulate organic P and particulate organic carbon can be excised without affecting our qualitative findings.From a starting point of eighteen nonlinear equations, we separately analysed a set of sixteen (approximately) linear equations which govern carbon and phosphorus dynamics and a pair of equations which explain the chemistry of the oxic deep ocean, the chemistry of the anoxic deep ocean and the nature of the transition between the two.
Having partitioned the system into two parts, we can elucidate the nature of the transition between oxic and anoxic oceans.A small change in system parameters produces abrupt, almost discontinuous, switches in the equilibrium concentrations of oxygen and reduced substances.We link this sensitivity in the model to the functional form prescribed for the removal of reduced substances as solid phases and the functional form for microbial respiration.Allison and Martiny (2008) refer to this kind of microbial model as a "black box" with "microorganisms buried within equation structure as kinetic constants and response functions".Our analysis highlights the need to compare the predictions of such studies with those of models that explicitly incorporate microbial biomass, in order to enhance understanding of how anoxia occurs.
With the nature of transition to anoxia established, we sought to determine the system parameters responsible for driving such a transition.Unfortunately, due to the scope of the original model (containing 69 parameters), the critical controlling parameter defies succinct characterisation.However, by focusing on a small subset of the system parameters (i.e., mixing rate and riverine input), we demonstrated that one 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 work could 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 variation (or covariation) in a few parameters of interest.More generally, we suggest that this article demonstrates the viability of adopting a systematic, mathematical approach in studying the behaviour of large biogeochemical models.The deduction of parameterised steady-state concentrations and equilibration time scales, as we have presented here, is generally beyond the reach of a purely computational approach to biogeochemistry.The dimensionless parameters of the resulting model (2.9) are given by the following three sets.First the O(1) parameters λ i : Next, the small, but not very small, parameters δ i : Finally, the very small parameters ε i : When the rescaling introduced before (2.16) is done, the new dimensionless parameters are given by

Dimensional parameters
In this section, we present the definitions and values of all dimensional constants associated with the model (2.1).Before we do so, we must first document the parameters of the original Slomp model.The values assigned to the parameters of the Slomp model, as well as their units and physical interpretations, are listed in tables 2 and 3.    Thus, our efforts to write an explicit formula for λ 6 s approx 3 − ν lead to extremely complicated formulae having no apparent simplification; it thus appears that the simple controlling parameters of the solution depend in a very complicated way on many of the physically prescribed parameters, and this dependence needs to be elucidated computationally (see Online Resource 1).

Figure 2 :
Figure2: Equilibration of soluble reactive phosphorus concentrations in each oceanic box.All variables are presented in dimensional form and time has been logarithmically 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.The initial value of each of the eighteen system variables was set to be 1 and (2.16) was solved numerically.

Figure 3 :
Figure3: The approximate slow manifold g = G(r), or g nullcline, (3.12), using parameter values from the Slomp model and ν set to 1 to represent a poorly-mixed ocean.In the lower curve, s 3 = 0.782 whereas in the upper curve, s 3 = 0.1.The variables g and r are dimensionless.

Figure 4 :
Figure4: ε 37 ṙ as a function of r given by (3.13) using parameter values from the Slomp model and ν set to 1 to represent a poorly-mixed ocean.In the upper curve, s 3 = 0.782 whereas in the lower curve, s 3 = 0.1.The quantities ε 37 ṙ and r are dimensionless.

Figure 5 :
Figure5: Steady-state values of r and g as a function of the dimensionless oceanic mixing parameter ν.Solutions of (2.16) were obtained numerically before variables were converted to their dimensional forms.

Figure 6
Figure 6: ṙ as a function of r plotted at ν = 4.1457 (blue curve) and ν = 4.1496 (black curve).Plotted variables are in dimensional form.The marked points correspond to steady-state values of r and all other system variables are at their numerically obtained steady-state values.

Figure 8 :
Figure8: Anoxia-onset metrics as a function of the mixing parameter ν.The quantity λ 6 s 3 −ν is presented in two forms, one using a numerical estimate of the s 3 steady-state value (λ 6 s num 3 − ν) and one using an analytical approximation (λ 6 s approx

Figure 9 :
Figure 9: Anoxia onset metric as a function of the riverine input parameter a 53 where we have set ν = 5.Positive values of λ 6 s 3 − ν correspond to anoxic conditions and negative values to oxic.

Table 1 :
F 5 ) and coastal upwelling (W F 6 ) are defined empirically, via constants that we will refer to as W k 1 , W k 5 and W k 6 respectively.Changes in circulation are modelled by multiplying the oceanic and coastal upwelling constants by the nondimensional parameters v o and v Definition of water fluxes in the Slomp model.The values of the constants W k 1 , W k 5 , W k 6 , v o and v c are given in table 2 of the appendix.2TheSlompandvan Capellen modelThe Slomp model divides the ocean into four distinct boxes: proximal coastal, distal coastal, surface ocean and deep ocean, having volumes W 1 -W 4 .Volume fluxes between the boxes are denoted by W F i , i = 1, 2, ..., 7. The boxes and fluxes are shown in figure1.As shown in table 1, the fluxes corresponding to river input (W F 1 ), ocean upwelling (W c , respectively.The remaining four fluxes in the oceanic circulation system then arise by imposing conservation of stationary water volume within each box.The values assigned to the water-cycle parameters are listed in table 2 of the appendix.

Table 2 :
Definition of the parameters of the Slomp model (continued below)