using AIBECS
using PyPlot, PyCall
using LinearAlgebra
using GR, Interact
The ideal tool for exploring global marine biogeochemical cycles
Algebraic Implicit Biogeochemical Elemental Cycling System
Check it on GitHub (look for AIBECS.jl)
Features:
Features (at present)
- no GUI but simple to use and fast
- nonlinear systems
- multiple tracers
- swappable circulations (not just the OCIM)
- made for parameter estimation/optimization and sensitivity analysis
To build a BGC model with the AIBECS, you just need to
1. Specify the local sources and sinks
2. Chose the ocean circulation
(3. Specify the particle transport, if any)
The 3D ocean grid is rearranged
into a 1D column vector.
And linear operators are represented by matrices, allowing use of powerful linear algebra tools.
Dissolved inorganic phosphrous (DIP)
(transported by the ocean circulation)
and particulate organic phosphrous (POP)
(transported by sinking particles)
To use AIBECS, we must recast each tracer equation for $\boldsymbol{x}_k$ into the generic form:
$$\frac{\partial \boldsymbol{x}_k}{\partial t} + \color{RoyalBlue}{\underbrace{\mathbf{T}(\boldsymbol{p})}_\textrm{transport}} \, \boldsymbol{x} = \color{ForestGreen}{\underbrace{\boldsymbol{G}(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N, \boldsymbol{p})}_\textrm{local sources and sinks}}$$where $\boldsymbol{p} =$ vector of model parameters
For DIP
$$\left[\frac{\partial}{\partial t} + \color{RoyalBlue}{\nabla \cdot (\boldsymbol{u} + \mathbf{K}\cdot\nabla )}\right] x_\mathsf{DIP} = \color{forestgreen}{-U(x_\mathsf{DIP}) + R(x_\mathsf{POP})},$$becomes
$$\left[\frac{\partial}{\partial t} + \color{RoyalBlue}{\mathbf{T}_{\mathsf{DIP}}(\boldsymbol{p})}\right] \boldsymbol{x}_\mathsf{DIP} = \color{forestgreen}{G_\mathsf{DIP}(\boldsymbol{x}_\mathsf{DIP}, \boldsymbol{x}_\mathsf{POP}, \boldsymbol{p})}$$For POP
$$\left[\frac{\partial}{\partial t} + \color{DarkOrchid}{\nabla \cdot \boldsymbol{w}}\right] x_\mathsf{POP} = \color{forestgreen}{U(x_\mathsf{DIP}) - R(x_\mathsf{POP})}.$$becomes
$$\left[\frac{\partial}{\partial t} + \color{RoyalBlue}{\mathbf{T}_{\mathsf{POP}}(\boldsymbol{p})}\right] \boldsymbol{x}_\mathsf{POP} = \color{forestgreen}{G_\mathsf{POP}(\boldsymbol{x}_\mathsf{DIP}, \boldsymbol{x}_\mathsf{POP}, \boldsymbol{p})}$$For DIP, the advective–eddy-diffusive transport operator, $\nabla \cdot (\boldsymbol{u} + \mathbf{K}\cdot\nabla)$, is converted into the OCIM1 transport matrix T
, which can be loaded via the OCIM1.load
function:
wet3D, grd, T = OCIM1.load();
We assign T
as the transport for DIP this way
T_DIP(p) = T
The matrix $\mathbf{T}$ acts on the column vector of the tracers.
What does T
look like?
T
For POP, the particle flux divergence (PFD) operator, $\nabla \cdot \boldsymbol{w}$, is created via buildPFD
. We assign the transport for POP via
T_POP(p) = buildPFD(grd, wet3D, sinking_speed = w(p))
The settling velocity, w(p)
, is assumed linearly increasing with depth z
to yield a "Martin curve profile"
w(p) = p.w₀ .+ p.w′ * (z .|> ustrip)
The vector of depths, z
, is defined through iwet
, which converts from 3D to 1D.
iwet = findall(wet3D)
z = grd.depth_3D[iwet]
relu(x) = (x .≥ 0) .* x
zₑ = 80u"m" # depth of the euphotic zone
function uptake(DIP, p)
τ, k = p.τ, p.k
DIP⁺ = relu(DIP)
return @. 1/τ * DIP⁺^2 / (DIP⁺ + k) * (z ≤ zₑ)
end
function plot_uptake()
PyPlot.figure(figsize=(8,4))
k = 0.5
τ = 5.0
x = 0:0.01:1.5
y = [1/τ*x^2/(x+k) for x in x]
y2 = [1/τ*x for x in x]
PyPlot.plot(x,y)
PyPlot.plot(x,y2)
PyPlot.plot(k*one.(x),y, ":")
PyPlot.xlabel("DIP (mol m⁻³)")
PyPlot.ylabel("uptake")
PyPlot.legend(("uptake (mol m⁻³ s⁻¹)", "1/τ DIP", "k (=$k mol m⁻³)"))
PyPlot.title("Uptake rate")
end
plot_uptake()
remineralization(POP, p) = p.κ * POP
geores(x, p) = (p.xgeo .- x) / p.τgeo
G_DIP(DIP, POP, p) = -uptake(DIP, p) + remineralization(POP, p) + geores(DIP, p)
G_POP(DIP, POP, p) = uptake(DIP, p) - remineralization(POP, p)
First, create a table of parameters using the AIBECS interface
t = empty_parameter_table() # empty table of parameters
add_parameter!(t, :xgeo, 2.12u"mmol/m^3")
add_parameter!(t, :τgeo, 1.0u"Myr")
add_parameter!(t, :k, 6.62u"μmol/m^3")
add_parameter!(t, :w₀, 0.64u"m/d")
add_parameter!(t, :w′, 0.13u"1/d")
add_parameter!(t, :κ, 0.19u"1/d")
add_parameter!(t, :τ, 236.52u"d")
t
We then create the parameters vector p
via
initialize_Parameters_type(t, "Pcycle_Parameters") # Generate the parameter type
p = Pcycle_Parameters()
The "state function" $\boldsymbol{F}$ is the total rate of change of the whole system (i.e., both DIP and POP) due to local sources and sinks and transport:
$$\frac{\partial \boldsymbol{x}}{\partial t} = \color{Brown}{\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})}$$For this coupled DIP–POP system,
$$\boldsymbol{F}\left(\begin{bmatrix}\boldsymbol{x}_\mathsf{DIP} \\ \boldsymbol{x}_\mathsf{POP}\end{bmatrix}, \boldsymbol{p}\right) = \begin{bmatrix}\color{forestgreen}{\boldsymbol{G}_\mathsf{DIP}(\boldsymbol{x}_\mathsf{DIP}, \boldsymbol{x}_\mathsf{POP}, \boldsymbol{p})} - \color{royalblue}{\boldsymbol{T}_\mathsf{DIP}} \, \boldsymbol{x}_\mathsf{DIP} \\ \color{forestgreen}{\boldsymbol{G}_\mathsf{POP}(\boldsymbol{x}_\mathsf{DIP}, \boldsymbol{x}_\mathsf{POP}, \boldsymbol{p})} - \color{royalblue}{\boldsymbol{T}_\mathsf{POP}} \, \boldsymbol{x}_\mathsf{POP} \end{bmatrix}$$We generate F
and ∇ₓF
via
nb = length(iwet)
F, ∇ₓF = state_function_and_Jacobian((T_DIP,T_POP), (G_DIP,G_POP), nb) ;
∇ₓF
?¶The Jacobian matrix is $\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x},\boldsymbol{p}) = \left[\frac{\partial F_i}{\partial x_j}\right]_{i,j}$, is useful for
With AIBECS, the Jacobian is computed automatically for you under the hood... using dual numbers!
We can directly solve for the steady-state, $\boldsymbol{s}$,
(using CTKAlg()
, a quasi-Newton root-finding algorithm from C. T. Kelley)
i.e., such that $\boldsymbol{F}(\boldsymbol{s},\boldsymbol{p}) = 0$:
x = ones(2nb)
prob = SteadyStateProblem(F, ∇ₓF, x, p)
s = solve(prob, CTKAlg(), preprint=" ").u
Plot DIP
DIP, POP = state_to_tracers(s, nb, 2)
DIP_unit = u"mmol/m^3"
lon, lat = grd.lon |> ustrip, grd.lat |> ustrip
function horizontal_slice(x, levels, title)
x_3D = fill(NaN, size(grd))
x_3D[iwet] .= x .|> ustrip
GR.figure(figsize=(10,5))
@manipulate for iz in 1:size(grd)[3]
GR.xlim([0,360])
GR.title(string(title, " at $(AIBECS.round(grd.depth[iz])) depth"))
GR.contourf(lon, lat, x_3D[:,:,iz]', levels=levels)
end
end
horizontal_slice(DIP * u"mol/m^3" .|> DIP_unit, 0:0.3:3, "P-cycle model: DIP [mmol m^{-3}]")
Plot POP
function zonal_slice(x, levels, title)
x_3D = fill(NaN, size(grd))
x_3D[iwet] .= x .|> ustrip
GR.figure(figsize=(10,5))
@manipulate for longitude in (grd.lon .|> ustrip)
GR.title(string(title, " at $(round(longitude))°"))
ilon = findfirst(ustrip.(grd.lon) .== longitude)
GR.contourf(lat, reverse(-grd.depth .|> ustrip), reverse(x_3D[:,ilon,:], dims=2), levels=levels)
end
end
zonal_slice(POP .|> relu .|> log10, -10:1:10, "P-cycle model: POP [log\\(mol m^{-3}\\)]")
We can also make publication-quality plots using Python's Cartopy
ccrs = pyimport("cartopy.crs")
cfeature = pyimport("cartopy.feature")
lon_cyc = [lon; 360+lon[1]]
DIP_3D = fill(NaN, size(grd)); DIP_3D[iwet] = DIP * u"mol/m^3" .|> u"mmol/m^3" .|> ustrip
DIP_3D_cyc = cat(DIP_3D, DIP_3D[:,[1],:], dims=2)
f1 = PyPlot.figure(figsize=(12,5))
@manipulate for iz in 1:24
withfig(f1, clear=true) do
ax = PyPlot.subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))
ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#CCCCCC") # gray land
plt = PyPlot.contourf(lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=0:0.1:3.5,
transform=ccrs.PlateCarree(), zorder=-1, extend="both")
plt2 = PyPlot.contour(plt, lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=plt.levels[1:5:end],
transform=ccrs.PlateCarree(), zorder=0, extend="both", colors="k")
ax.clabel(plt2, fmt="%2.1f", colors="k", fontsize=14)
cbar = PyPlot.colorbar(plt, orientation="vertical", extend="both", ticks=plt2.levels)
cbar.add_lines(plt2)
cbar.ax.tick_params(labelsize=14)
cbar.set_label(label="mmol / m³", fontsize=16)
PyPlot.title("Publication-quality DIP with Cartopy; depth = $(string(round(typeof(1u"m"),grd.depth[iz])))", fontsize=20)
end
end