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, z -> w(z,p))

T_POP (generic function with 1 method)

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

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

w (generic function with 1 method)

Uptake (DIP → POP)

For the uptake, $U$, we write

z = depthvec(grd)
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 DIP_geo, τ_geo = p
    return @. -$U(DIP,p) + $R(POP,p) + (DIP_geo - DIP) / τ_geo
end
function G_POP(DIP, POP, p)
    @unpack τ_geo = p
    return @. $U(DIP,p) - $R(POP,p) - POP / τ_geo
end

G_POP (generic function with 1 method)

where we have imposed a slow restoring of DIP to the global mean DIP_geo 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
using Unitful: m, d, s, yr, Myr, mol, mmol, μmol, μM
@initial_value @units struct PmodelParameters{U} <: AbstractParameters{U}
    w₀::U       |  0.64 | m/d
    w′::U       |  0.13 | m/d/m
    τ_DIP::U    | 230.0 | d
    k::U        |  6.62 | μmol/m^3
    z₀::U       |  80.0 | m
    τ_POP::U    |   5.0 | d
    τ_geo::U    |   1.0 | Myr
    DIP_geo::U  |  2.12 | mmol/m^3
end

initial_value (generic function with 24 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        FreeUnit…
─────┼────────────────────────────────────────────
   1 │ w₀          0.64           0.64  m d⁻¹
   2 │ w′          0.13           0.13  d⁻¹
   3 │ τ_DIP     230.0          230.0   d
   4 │ k           6.62           6.62  μmol m⁻³
   5 │ z₀         80.0           80.0   m
   6 │ τ_POP       5.0            5.0   d
   7 │ τ_geo       1.0            1.0   Myr
   8 │ DIP_geo     2.12           2.12  mmol m⁻³

We generate the state function F,

nb = sum(iswet(grd))
F = AIBECSFunction((T_DIP, T_POP), (G_DIP, G_POP), nb)

(::SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, AIBECS.var"#f#58"{Tuple{typeof(Main.T_DIP), typeof(Main.T_POP)}, UnitRange{Int64}, AIBECS.var"#G#56"{Tuple{typeof(Main.G_DIP), typeof(Main.G_POP)}, AIBECS.var"#tracers#54"{Int64, Int64}}, AIBECS.var"#tracer#55"{Int64, Int64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, AIBECS.var"#jac#63"{AIBECS.var"#T#60"{Tuple{typeof(Main.T_DIP), typeof(Main.T_POP)}, Int64, UnitRange{Int64}}, AIBECS.var"#∇ₓG#59"{Tuple{typeof(Main.G_DIP), typeof(Main.G_POP)}, Int64, Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}) (generic function with 1 method)

generate the steady-state problem,

@unpack DIP_geo = p
x = DIP_geo * ones(2nb) # initial guess
prob = SteadyStateProblem(F, x, p)

SteadyStateProblem with uType Vector{Float64}. In-place: false
u0: 382338-element Vector{Float64}:
 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

sol = solve(prob, CTKAlg()).u

382338-element Vector{Float64}:
 0.0020967060231664295
 0.0021333454639532027
 0.0018399997020561354
 0.0017256341824069535
 0.001612328585777842
 0.0015535540266685076
 0.001458541565317635
 0.001570008107294452
 0.001457342544838417
 0.00135850268497681
 ⋮
 2.7883468616763368e-9
 1.6701226349887806e-9
 1.8293710129060015e-9
 2.873136469222496e-9
 2.993647120876417e-9
 2.8588624377785625e-9
 3.2379918118282254e-9
 3.0900967665459405e-9
 3.1845337501702697e-9

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

DIP, POP = state_to_tracers(sol, grd) # unpack tracers

([0.0020967060231664295, 0.0021333454639532027, 0.0018399997020561354, 0.0017256341824069535, 0.001612328585777842, 0.0015535540266685076, 0.001458541565317635, 0.001570008107294452, 0.001457342544838417, 0.00135850268497681  …  0.001411275166907556, 0.001412264306571928, 0.0013911523315427185, 0.0013911433495659214, 0.0021377405632125503, 0.0021393993335517133, 0.002135737129869365, 0.0021249929534287804, 0.0021344756362701237, 0.0021379362067694492], [2.6134960067708514e-5, 2.6593099800802548e-5, 2.292510803375331e-5, 2.1495086727744334e-5, 2.007832138108026e-5, 1.934340962415626e-5, 1.8155384587571548e-5, 1.9549149860286034e-5, 1.8140392190430782e-5, 1.6904513204783926e-5  …  3.0408296778339733e-9, 2.7883468616763368e-9, 1.6701226349887806e-9, 1.8293710129060015e-9, 2.873136469222496e-9, 2.993647120876417e-9, 2.8588624377785625e-9, 3.2379918118282254e-9, 3.0900967665459405e-9, 3.1845337501702697e-9])

First, let's look at the mean profile

using Plots
plothorizontalmean(DIP * (mol/m^3) .|> μM, grd)

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

plothorizontalslice(DIP * (mol/m^3) .|> μM, grd, depth=1000m, color=:viridis)

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

plotverticalintegral(U(DIP,p) * (mol/m^3/s) .|> mmol/yr/m^3, grd, color=:algae)

Or look at what is exported below 500 m

plothorizontalslice(POP .* w(z,p) * (mol/m^3*m/s) .|> mmol/yr/m^2, grd, depth=500m, color=:inferno, rev=true)

Now let's make our model a little fancier and use a fine topographic map to refine the remineralization profile. For this, we will use the ETOPO dataset, which can be downloaded by AIBECS via

f_topo = ETOPO.fractiontopo(grd)

191169-element Vector{Float64}:
 0.0006887052341597796
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.9904166666666666
 0.9914583333333333
 0.9657638888888889
 0.9586805555555555
 0.9751388888888889
 0.9947222222222222
 0.9328472222222223
 0.9486111111111111
 0.9800694444444444

f_topo is the fraction of sediments located in each wet box of the grd grid. We can use it to redefine the transport operator for sinking particles to take into consideration the subgrid topography, such that the fine-resolution sediments intercept settling POP.

T_POP2(p) = transportoperator(grd, z -> w(z,p); frac_seafloor=f_topo)

T_POP2 (generic function with 1 method)

With this new vertical transport for POP, we can recreate our problem, solve it again

F2 = AIBECSFunction((T_DIP, T_POP2), (G_DIP, G_POP), nb)
prob2 = SteadyStateProblem(F2, x, p)
sol2 = solve(prob2, CTKAlg()).u
DIP2, POP2 = state_to_tracers(sol2, grd) # unpack tracers

([0.002098085071985876, 0.0021343221222590372, 0.001840983313044677, 0.001726670590169758, 0.0016137061456333046, 0.0015554917832502086, 0.0014607122946576803, 0.0015728691982204667, 0.0014603298577997833, 0.0013625583518128337  …  0.0014213554906768006, 0.0014223553341716668, 0.001401365142005268, 0.001401359910709824, 0.0021407469293816704, 0.002142403108241433, 0.0021387585449112196, 0.0021280771897256535, 0.002137500584097961, 0.0021409435430922575], [2.615985726077029e-5, 2.66053119427593e-5, 2.2937407071766355e-5, 2.150804591211601e-5, 2.009554627567894e-5, 1.9367639141425456e-5, 1.818252710802236e-5, 1.958492467288954e-5, 1.8177745166624564e-5, 1.6955224585933824e-5  …  3.1092058441255798e-9, 2.8716550526604888e-9, 1.8040447193459132e-9, 1.9772526509285164e-9, 2.937879236394193e-9, 3.0701301969675845e-9, 2.9462721482362614e-9, 3.2997240879145508e-9, 3.15904173541064e-9, 3.266028827607643e-9])

and check the difference

plotzonalaverage((DIP2 - DIP) ./ DIP .|> u"percent", grd, color=:balance, clim=(-5,5))

This zonal average shows how much DIP is redistributed as it is prevented from sinking out of the surface layers with the new subgrid parameterization.

Let's look at the vertical average.

plotverticalaverage((DIP2 - DIP) ./ DIP .|> u"percent", grd, color=:balance, clim=(-10,10))

This shows minor changes on the order of 1%, on the global scale, except along the coasts, which retain much more DIP with the subgrid topography.

This page was generated using Literate.jl.