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
┌ Warning: Over-writing registration of the datadep │ name = "AIBECS-OCIM0.1" └ @ DataDeps ~/.julia/packages/DataDeps/ae6dT/src/registration.jl:15 ┌ Info: You are about to use the OCIM0.1 model. │ If you use it for research, please cite: │ │ - Primeau, F. W., Holzer, M., and DeVries, T. (2013), Southern Ocean nutrient trapping and the efficiency of the biological pump, J. Geophys. Res. Oceans, 118, 2547–2564, doi:10.1002/jgrc.20181. │ - DeVries, T. and F. Primeau, 2011: Dynamically and Observationally Constrained Estimates of Water-Mass Distributions and Ages in the Global Ocean. J. Phys. Oceanogr., 41, 2381–2401, doi:10.1175/JPO-D-10-05011.1 │ │ │ You can find the corresponding BibTeX entries in the CITATION.bib file │ at the root of the AIBECS.jl package repository. └ (Look for the "DeVries_Primeau_2011" and "Primeau_etal_2013" keys.)
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)
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.)
For the remineralization, $R$, we write
function R(x,p)
@unpack τ_POP = p
return x / τ_POP
end
R (generic function with 1 method)
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.var"##354".T_DIP), typeof(Main.var"##354".T_POP)}, UnitRange{Int64}, AIBECS.var"#G#56"{Tuple{typeof(Main.var"##354".G_DIP), typeof(Main.var"##354".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.var"##354".T_DIP), typeof(Main.var"##354".T_POP)}, Int64, UnitRange{Int64}}, AIBECS.var"#∇ₓG#59"{Tuple{typeof(Main.var"##354".G_DIP), typeof(Main.var"##354".G_POP)}, Int64, Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), 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.0020967060231705586 0.0021333454639574046 0.0018399997020597568 0.0017256341824103496 0.0016123285857810134 0.0015535540266715638 0.0014585415653205036 0.0015700081072975408 0.0014573425448412827 0.0013585026849794804 ⋮ 2.7883468616818437e-9 1.6701226349920786e-9 1.8293710129096136e-9 2.8731364692281686e-9 2.993647120882326e-9 2.8588624377842093e-9 3.237991811834619e-9 3.0900967665520413e-9 3.1845337501765604e-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.0020967060231705586, 0.0021333454639574046, 0.0018399997020597568, 0.0017256341824103496, 0.0016123285857810134, 0.0015535540266715638, 0.0014585415653205036, 0.0015700081072975408, 0.0014573425448412827, 0.0013585026849794804 … 0.0014112751669103055, 0.0014122643065746792, 0.001391152331545428, 0.001391143349568631, 0.002137740563216725, 0.002139399333555891, 0.002135737129873535, 0.002124992953432931, 0.002134475636274292, 0.0021379362067736243], [2.6134960067760146e-5, 2.6593099800855085e-5, 2.292510803379859e-5, 2.14950867277868e-5, 2.007832138111992e-5, 1.9343409624194476e-5, 1.8155384587607418e-5, 1.954914986032466e-5, 1.8140392190466615e-5, 1.6904513204817316e-5 … 3.040829677839977e-9, 2.7883468616818437e-9, 1.6701226349920786e-9, 1.8293710129096136e-9, 2.8731364692281686e-9, 2.993647120882326e-9, 2.8588624377842093e-9, 3.237991811834619e-9, 3.0900967665520413e-9, 3.1845337501765604e-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)
┌ Warning: Binning ETOPO to grd. This may take a few seconds └ @ AIBECS.ETOPO ~/work/AIBECS.jl/AIBECS.jl/src/ETOPO.jl:120 ┌ Warning: Checksum not provided, add to the Datadep Registration the following hash line │ hash = "55ebef118c1bd4d8682464d4d24f0d0850d24160c996d1ff62ca496010ffdf61" └ @ DataDeps ~/.julia/packages/DataDeps/ae6dT/src/verification.jl:44 7-Zip (a) [64] 17.04 : Copyright (c) 1999-2021 Igor Pavlov : 2017-08-28 p7zip Version 17.04 (locale=utf8,Utf16=on,HugeFiles=on,64 bits,3 CPUs x64) Scanning the drive for archives: 1 file, 401881655 bytes (384 MiB) Extracting archive: /Users/runner/.julia/datadeps/ETOPO_bedrock/ETOPO1_Bed_g_gdal.grd.gz -- Path = /Users/runner/.julia/datadeps/ETOPO_bedrock/ETOPO1_Bed_g_gdal.grd.gz Type = gzip Headers Size = 32 Everything is Ok Size: 933250996 Compressed: 401881655 ┌ Info: You are about to use the ETOPO data set. │ If you use it for research, please cite: │ │ - Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M [access date]. │ │ You can find the corresponding BibTeX entries in the CITATION.bib file │ at the root of the AIBECS.jl package repository. └ (Look for the "Amante_Eakins_2009" key.)
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.002098085071990412, 0.002134322122263653, 0.0018409833130486558, 0.001726670590173489, 0.0016137061456367907, 0.0015554917832535692, 0.001460712294660835, 0.0015728691982238652, 0.0014603298578029374, 0.0013625583518157756 … 0.0014213554906798663, 0.0014223553341747345, 0.0014013651420082711, 0.0014013599107128268, 0.002140746929386296, 0.0021424031082460625, 0.002138758544915841, 0.0021280771897302523, 0.0021375005841025797, 0.002140943543096884], [2.6159857260827024e-5, 2.6605311942817016e-5, 2.2937407071816106e-5, 2.1508045912162667e-5, 2.0095546275722533e-5, 1.936763914146748e-5, 1.8182527108061812e-5, 1.958492467293204e-5, 1.817774516666401e-5, 1.6955224585970612e-5 … 3.1092058441323234e-9, 2.8716550526667183e-9, 1.8040447193498278e-9, 1.977252650932807e-9, 2.9378792364005663e-9, 3.070130196974248e-9, 2.9462721482426546e-9, 3.2997240879217117e-9, 3.1590417354174945e-9, 3.2660288276147305e-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 notebook was generated using Literate.jl.