A coupled PO₄–POP model

In this tutorial we will explicitly simulate 2 tracers whose distributions control and feed back on each other.

We consider a simple model for the cycling of phosphorus with 2 state variables consisting of phosphate (PO₄) AKA dissolved inorganic phosphorus (DIP) and particulate organic phosphorus (POP). The dissolved phases are transported by advection and diffusion whereas the particulate phase sinks rapidly down the water column without any appreciable transport by the circulation.

The governing equations that couple the 3D concentration fields of DIP and POP, denoted $x_\mathsf{DIP}$ and $x_\mathsf{POP}$, respectively, are:

\[\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \mathbf{K}\nabla )\right] x_\mathsf{DIP} = -U(x_\mathsf{DIP}) + R(x_\mathsf{POP}),\]

and

\[\left[\frac{\partial}{\partial t} + \nabla \cdot \boldsymbol{w}\right] x_\mathsf{POP} = U(x_\mathsf{DIP}) - R(x_\mathsf{POP}).\]

The $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ and $\nabla \cdot \boldsymbol{w}$ operators represent the ocean circulation and the sinking of particles, respectively. (Tracer transport operators are described in the documentation.)

The function $U$ represents the biological uptake of DIP by phytoplankton, which we model here as

\[U(x_\mathsf{DIP}) = \frac{x_\mathsf{DIP}}{\tau_\mathsf{DIP}} \, \frac{x_\mathsf{DIP}}{x_\mathsf{DIP} + k} \, (z < z_0),\]

with the timescale, $\tau$, the half-saturation rate $k$, and the depth $z_0$ as parameters.

The function $R$ defines the remineralization rate of POP, which converts POP back into DIP. For the remineralization, we simply use a linear rate constant, i.e.,

\[R(x_\mathsf{POP}) = \frac{x_\mathsf{POP}}{\tau_\mathsf{POP}}.\]

We start by telling Julia we want to use the AIBECS and the OCIM0.1 circulation for DIP.

using AIBECS
grd, T_OCIM = OCIM0.load()
T_DIP(p) = T_OCIM
T_DIP (generic function with 1 method)

For the sinking of particles, we use the transportoperator function

T_POP(p) = transportoperator(grd, w = w(p))
T_POP (generic function with 1 method)

for which we need to define the sinking speed w(p) as a function of the parameters p. Following the assumption that $w = w_0 + w' z$ increases linearly with depth, we write it as

function w(p)
    @unpack w₀, w′ = p
    return @. w₀ + w′ * z
end
w (generic function with 1 method)

For this to work, we must create a vector of depths, z, which is simply done via

z = depthvec(grd)
191169-element Array{Float64,1}:
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
    ⋮              
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
Uptake (DIP → POP)

For the uptake, $U$, we write

function U(x,p)
    @unpack τDIP, k, z₀ = p
    return @. x/τDIP * x/(x+k) * (z≤z₀) * (x≥0)
end
U (generic function with 1 method)

where we have "unpacked" the parameters to make the code clearer and as close to the mathematical equation as possible. (Note we have also added a constraint that x must be positive for uptake to happen.)

Remineralization (POP → DIP)

For the remineralization, $R$, we write

function R(x,p)
    @unpack τPOP = p
    return x / τPOP
end
R (generic function with 1 method)
Net sources and sinks

We lump the sources and sinks into G functions for DIP and POP.

function G_DIP(DIP, POP, p)
    @unpack D̅I̅P̅, τgeo = p
    return @. -$U(DIP,p) + $R(POP,p) + (D̅I̅P̅ - DIP) / τgeo
end
function G_POP(DIP, POP, p)
    return U(DIP,p) - R(POP,p)
end
G_POP (generic function with 1 method)

where we have imposed a slow restoring of DIP to the global mean D̅I̅P̅ to prescribe the global mean concentration. (The $ signs in front of U and R protect them from the broadcast macro @.)

We now define and build the parameters.

In this tutorial we will specify some initial values for the parameters and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@units @initial_value struct PmodelParameters{U} <: AbstractParameters{U}
    w₀::U   |  0.64 | u"m/d"
    w′::U   |  0.13 | u"m/d/m"
    τDIP::U | 230.0 | u"d"
    k::U    |  6.62 | u"μmol/m^3"
    z₀::U   |  80.0 | u"m"
    τPOP::U |   5.0 | u"d"
    τgeo::U |   1.0 | u"Myr"
    D̅I̅P̅::U  |  2.12 | u"mmol/m^3"
end
units (generic function with 47 methods)

Finally, thanks to the initial values we provided, we can instantiate the parameter vector succintly as

p = PmodelParameters()

│ Row │ Symbol │ Value   │ Initial value │ Unit      │
│     │ Symbol │ Float64 │ Float64       │ Unitful…  │
├─────┼────────┼─────────┼───────────────┼───────────┤
│ 1   │ w₀     │ 0.64    │ 0.64          │ m d^-1    │
│ 2   │ w′     │ 0.13    │ 0.13          │ d^-1      │
│ 3   │ τDIP   │ 230.0   │ 230.0         │ d         │
│ 4   │ k      │ 6.62    │ 6.62          │ μmol m^-3 │
│ 5   │ z₀     │ 80.0    │ 80.0          │ m         │
│ 6   │ τPOP   │ 5.0     │ 5.0           │ d         │
│ 7   │ τgeo   │ 1.0     │ 1.0           │ Myr       │
│ 8   │ D̅I̅P̅    │ 2.12    │ 2.12          │ mmol m^-3 │

We generate the state function F and its Jacobian ∇ₓF,

nb = sum(iswet(grd))
F, ∇ₓF = state_function_and_Jacobian((T_DIP, T_POP), (G_DIP, G_POP), nb)
(AIBECS.var"#F#27"{Tuple{typeof(Main.ex-3_Pmodel.T_DIP),typeof(Main.ex-3_Pmodel.T_POP)},AIBECS.var"#tracer#22"{Int64,Int64},AIBECS.var"#G#25"{Tuple{typeof(Main.ex-3_Pmodel.G_DIP),typeof(Main.ex-3_Pmodel.G_POP)},AIBECS.var"#tracers#21"{Int64,Int64}}}((Main.ex-3_Pmodel.T_DIP, Main.ex-3_Pmodel.T_POP), AIBECS.var"#tracer#22"{Int64,Int64}(191169, 2), AIBECS.var"#G#25"{Tuple{typeof(Main.ex-3_Pmodel.G_DIP),typeof(Main.ex-3_Pmodel.G_POP)},AIBECS.var"#tracers#21"{Int64,Int64}}((Main.ex-3_Pmodel.G_DIP, Main.ex-3_Pmodel.G_POP), AIBECS.var"#tracers#21"{Int64,Int64}(191169, 2))), AIBECS.var"#∇ₓF#30"{AIBECS.var"#T#23"{Tuple{typeof(Main.ex-3_Pmodel.T_DIP),typeof(Main.ex-3_Pmodel.T_POP)}},AIBECS.var"#∇ₓG#29"{Tuple{typeof(Main.ex-3_Pmodel.G_DIP),typeof(Main.ex-3_Pmodel.G_POP)},Int64,Int64}}(AIBECS.var"#T#23"{Tuple{typeof(Main.ex-3_Pmodel.T_DIP),typeof(Main.ex-3_Pmodel.T_POP)}}((Main.ex-3_Pmodel.T_DIP, Main.ex-3_Pmodel.T_POP)), AIBECS.var"#∇ₓG#29"{Tuple{typeof(Main.ex-3_Pmodel.G_DIP),typeof(Main.ex-3_Pmodel.G_POP)},Int64,Int64}((Main.ex-3_Pmodel.G_DIP, Main.ex-3_Pmodel.G_POP), 191169, 2)))

generate the steady-state problem,

@unpack D̅I̅P̅ = p
x = D̅I̅P̅ * ones(2nb) # initial guess
prob = SteadyStateProblem(F, ∇ₓF, x, p)
SteadyStateProblem with uType Array{Float64,1}
u0: [0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004  …  0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004, 0.0021200000000000004]

and solve it

s = solve(prob, CTKAlg()).u
382338-element Array{Float64,1}:
 0.002056497764443025 
 0.0020749147922378467
 0.001729268855455424 
 0.0016232763386096164
 0.001489871068305587 
 0.0014351512033222947
 0.0013340342196832852
 0.0014636961680406385
 0.001327345002371588 
 0.001226108004021758 
 ⋮                    
 3.739881068755677e-9 
 2.1533754283537364e-9
 2.3168285890004067e-9
 3.906941812936834e-9 
 4.115361931140372e-9 
 3.938340795598846e-9 
 4.42333901595675e-9  
 4.262506794662443e-9 
 4.398742297227606e-9 

We can look at different the DIP and POP fields using the Plots.jl recipes.

DIP, POP = state_to_tracers(s, grd) # unpack tracers
([0.002056497764443025, 0.0020749147922378467, 0.001729268855455424, 0.0016232763386096164, 0.001489871068305587, 0.0014351512033222947, 0.0013340342196832852, 0.0014636961680406385, 0.001327345002371588, 0.001226108004021758  …  0.0013577912362383716, 0.0013591449170113992, 0.0013287055936139585, 0.0013286920404612425, 0.0022054153729531875, 0.002207445941740303, 0.0022030371074698452, 0.002190172786064569, 0.002201506111467498, 0.0022056665407433123], [2.1489486144294227e-5, 2.168255342399456e-5, 1.8059124436327874e-5, 1.6948001676310133e-5, 1.5549513657990965e-5, 1.4975886424524809e-5, 1.3915881641085637e-5, 1.5275122601822942e-5, 1.3845759005895051e-5, 1.2784500335188343e-5  …  4.082888958842055e-9, 3.739881068755677e-9, 2.1533754283537364e-9, 2.3168285890004067e-9, 3.906941812936834e-9, 4.115361931140372e-9, 3.938340795598846e-9, 4.42333901595675e-9, 4.262506794662443e-9, 4.398742297227606e-9])

We can plot the concentration of DIP at a given depth via, e.g.,

using Plots
horizontalslice(DIP * u"mol/m^3" .|> u"μM", grd, depth=1000u"m", color=:viridis)
50 100 150 200 250 300 350 -60 -30 0 30 60 Longitude (°) Latitude (°) 0.5 1.0 1.5 2.0 2.5 3.0 3.5 ?M

Or have a look at a map of the uptake at the surface

verticalintegral(U(DIP,p) * u"mol/m^3/s" .|> u"mmol/yr/m^3", grd, color=:algae)
50 100 150 200 250 300 350 -60 -30 0 30 60 Longitude (°) Latitude (°)