{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "using AIBECS\n", "using PyPlot, PyCall\n", "using LinearAlgebra\n", "using GR, Interact" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "

AIBECS.jl

\n", "\n", "*The ideal tool for exploring global marine biogeochemical cycles*\n", "\n", "**A**lgebraic **I**mplicit **B**iogeochemical **E**lemental **C**ycling **S**ystem\n", "\n", "Check it on GitHub (look for [AIBECS.jl](https://github.com/briochemc/AIBECS.jl))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "\n", "\n", "
A Julia package developed by Benoît Pasquier at the Department of Earth System Sciences, UCI
with François Primeau and J. Keith Moore
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Outline\n", "\n", "1. Motivation and concept\n", "1. Example 1: Radiocarbon\n", " 1. Toy model circulation\n", " 1. OCIM1\n", "1. Example 2: Phosphorus cycle\n", "1. AIBECS ecosystem" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Motivation: Starting from the AWESOME OCIM\n", "\n", "\n", "\n", "The AWESOME OCIM (for A Working Environment for Simulating Ocean Movement and Elemental cycling in an Ocean Circulation Inverse Model framework) by **Seth John** (USC)\n", "\n", "Features: **GUI**, **simple to use**, **fast**, and **good circulation**
(thanks to the OCIM1 by **Tim DeVries** (UCSB))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Motivation: comes the AIBECS\n", "\n", "\n", "\n", "Features (at present)\n", "\n", "\\- **simple to use**
\n", "\\- **fast**
\n", "\\- **Julia** instead of MATLAB (free, open-source, better performance, and better syntax)
\n", "\\- **nonlinear** systems
\n", "\\- **multiple tracers**
\n", "\\- **Other circulations** (not just the OCIM)
\n", "\\- **Parameter estimation/optimization** and **Sensitivity analysis** (shameless plug: F-1 algorithm seminar tomorrow at the School of Mathematics)
\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "# AIBECS Concept: a simple interface\n", "\n", "To build a BGC model with the AIBECS, you just need to\n", "\n", "**1.** Specify the **local sources and sinks**\n", "\n", "**2.** Chose the **ocean circulation**
\n", "\n", "(**3.** Specify the **particle transport**, if any)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "# AIBECS concept: Vectorization\n", "\n", "The **3D ocean grid** is rearranged
into a **1D column vector**.\n", "\n", "And **linear operators** are represented by **matrices**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "

Example 1: Radiocarbon, a tracer for water age

\n", "\n", "
\n", "\n", "\n", "\n", "*Image credit: Luke Skinner, University of Cambridge*\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "

Tracer equation: transport + sources and sinks

\n", "\n", "The **Tracer equation**\n", " ($x=$ Radiocarbon concentration)\n", "\n", "$$\\frac{\\partial x}{\\partial t} + \\color{RoyalBlue}{\\nabla \\cdot \\left[ \\boldsymbol{u} - \\mathbf{K} \\cdot \\nabla \\right]} x = \\color{ForestGreen}{\\underbrace{\\Lambda(x)}_{\\textrm{air–sea exchange}} - \\underbrace{x / \\tau}_{\\textrm{radioactive decay}}}$$\n", "\n", "becomes \n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} + \\color{RoyalBlue}{\\mathbf{T}} \\, \\boldsymbol{x} = \\color{ForestGreen}{\\mathbf{\\Lambda}(\\boldsymbol{x}) - \\boldsymbol{x} / \\tau}.$$\n", "\n", "with the **transport matrix**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Translating to AIBECS Code is easy\n", "\n", "To use AIBECS, we must recast each tracer equation,\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} + \\color{RoyalBlue}{\\mathbf{T}} \\, \\boldsymbol{x} = \\color{ForestGreen}{\\mathbf{\\Lambda}(\\boldsymbol{x}) - \\boldsymbol{x} / \\tau}$$\n", "\n", "here, into the generic form:\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} + \\color{RoyalBlue}{\\mathbf{T}(\\boldsymbol{p})} \\, \\boldsymbol{x} = \\color{ForestGreen}{\\boldsymbol{G}(\\boldsymbol{x}, \\boldsymbol{p})}$$\n", "\n", "where $\\boldsymbol{p} =$ vector of model parameters\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Circulation 1: The 2×2×2 *Primeau* model\n", "\n", "\n", "\n", "\\- **ACC**: Antarctic Circumpolar Current\n", "\n", "\\- **MOC**: Meridional Overturning Circulation\n", "\n", "\\- High-latitude mixing\n", "\n", "\n", "(Credit: François Primeau, and Louis Primeau for the image)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Load the **circulation** via `load`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "wet3D, grd, T = Primeau_2x2x2.load();" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "`wet3D` is the mask of \"wet\" points" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "wet3D" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "`grd` is the grid of the circulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "grd" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "We can check the depth of the boxes arranged in 3D" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "grd.depth_3D" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "### The matrix $\\mathbf{T}$ acts on the column vector\n", "\n", "What does `T` look like?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "T" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "A sparse matrix is indexed by its non-zero values,
\n", "but we can check it out in full using `Matrix`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "Matrix(T)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Sources and sinks\n", "\n", "Tracer equation reminder:\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} + \\mathbf{T}(\\boldsymbol{p}) \\, \\boldsymbol{x} = \\boldsymbol{G}(\\boldsymbol{x}, \\boldsymbol{p})$$\n", "\n", "Let's write $\\boldsymbol{G}(\\boldsymbol{x}, \\boldsymbol{p}) = \\mathbf{\\Lambda}(\\boldsymbol{x}) - \\boldsymbol{x} / \\tau$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "G(x,p) = Λ(x,p) - x / p.τ" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Air–sea gas exchange\n", "\n", "And define the air–sea gas exchange $\\mathbf{\\Lambda}(\\boldsymbol{x}) = \\frac{\\lambda}{h} (R_\\mathsf{atm} - \\boldsymbol{x})$ at the surface with a piston velocity $\\lambda$ over the top layer of height $h$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "function Λ(x,p)\n", " λ, h, Ratm = p.λ, p.h, p.Ratm\n", " return @. λ / h * (Ratm - x) * (z == z₀)\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "Define `z` the depths in vector form.
\n", "(`iwet` converts from 3D to 1D)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "iwet = findall(wet3D)\n", "z = grd.depth_3D[iwet] " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Define `z₀` the depth of the top layer" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "z₀ = z[1]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "So that `z .== z₀` is `true` at the surface layer" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "z .== z₀" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Model parameters\n", "\n", "First, create a table of parameters using the AIBECS API" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "t = empty_parameter_table()\n", "add_parameter!(t, :τ, 5730u\"yr\" / log(2)) # radioactive decay e-folding timescale\n", "add_parameter!(t, :λ, 50u\"m\" / 10u\"yr\") # piston velocity\n", "add_parameter!(t, :h, grd.δdepth[1]) # top layer height\n", "add_parameter!(t, :Ratm, 1.0u\"mol/m^3\") # atmospheric concentration\n", "t " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Then, chose a name for the parameters (here `C14_parameters`), and create the vector `p`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "initialize_Parameters_type(t, \"C14_parameters\")\n", "p = C14_parameters() " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Note `p` has units! \n", "\n", "In AIBECS, you give your parameters units and they are **automatically converted to SI units** under the hood.\n", "\n", "(And they are converted back for pretty printing!)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# State function (and Jacobian)\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} = \\boldsymbol{G}(\\boldsymbol{x}, \\boldsymbol{p}) - \\mathbf{T}(\\boldsymbol{p}) \\, \\boldsymbol{x} = \\color{Brown}{\\boldsymbol{F}(\\boldsymbol{x}, \\boldsymbol{p})}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "We generate `F` and `∇ₓF` via" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "F, ∇ₓF = state_function_and_Jacobian(p -> T, G) ; " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## The state function `F(x,p)`\n", "\n", "Let's try `F` on a random state vector `x`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "x = 0.5p.Ratm * ones(5)\n", "F(x,p)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## The Jacobian `∇ₓF`\n", "\n", "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 \n", "- **implicit** time-steps\n", "- **solving** the **steady-state** system\n", "- **optimization** / **uncertainty analysis**\n", "\n", "With AIBECS, the **Jacobian** is **computed automatically** for you under the hood... using **dual numbers**!
\n", "(Come to my Applied seminar tomorrow for more on dual numbers and... hyperdual numbers!)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's try `∇ₓF` at `x`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "Matrix(∇ₓF(x,p))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Time stepping\n", "\n", "AIBECS provides schemes for time-stepping\n", "- Euler forward\n", "- Euler backward\n", "- **Crank-Nicolson**\n", "- Crank-Nicolson leap-frog" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's track the evolution of `x` through time\n", "\n", "Define a function to apply the time steps `n` times for a time span of `Δt` starting from `x₀`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "function time_steps(x₀, Δt, n, F, ∇ₓF)\n", " x_hist = [x₀]\n", " δt = Δt / n\n", " for i in 1:n\n", " push!(x_hist, AIBECS.crank_nicolson_step(last(x_hist), p, δt, F, ∇ₓF))\n", " end\n", " return reduce(hcat, x_hist), 0:δt:Δt\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's run the model for 5000 years starting with `x = 1` everywhere:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "Δt = 5000u\"yr\" |> u\"s\" |> ustrip\n", "x₀ = p.Ratm * ones(5) \n", "x_hist, t_hist = time_steps(x₀, Δt, 1000, F, ∇ₓF) " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "#### Plotting the output is easy\n", "\n", "The radiocarbon age, `C14age`, is given by $\\log(R_{\\mathrm{atm}}/\\boldsymbol{x}) \\tau$ because $\\boldsymbol{x}\\sim R_{\\mathrm{atm}} \\exp(-t/\\tau)$\n", "\n", "Let's plot its evolution with time:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "C14age_hist = log.(p.Ratm ./ x_hist) * (p.τ * u\"s\" |> u\"yr\" |> ustrip)\n", "PyPlot.figure(figsize=(8,4))\n", "PyPlot.plot(t_hist .* 1u\"s\" .|> u\"yr\" .|> ustrip, C14age_hist')\n", "PyPlot.xlabel(\"simulation time (years)\")\n", "PyPlot.ylabel(\"¹⁴C age (years)\")\n", "PyPlot.legend(\"box \" .* string.(findall(vec(wet3D))))\n", "PyPlot.title(\"Simulation of the evolution of ¹⁴C age with Crank-Nicolson time steps\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Solve directly for the steady state\n", "\n", "Instead, we can directly solve for the **steady-state**, $\\boldsymbol{s}$,
\n", "(using `CTKAlg()`, a quasi-Newton root-finding algorithm from C. T. Kelley)\n", "\n", "i.e., such that $\\boldsymbol{F}(\\boldsymbol{s},\\boldsymbol{p}) = 0$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "prob = SteadyStateProblem(F, ∇ₓF, x, p)\n", "s = solve(prob, CTKAlg()).u" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "gives the age" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "log.(p.Ratm ./ s) * (p.τ * u\"s\" |> u\"yr\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "# 35'000 years without the steady-state solver!\n", "\n", "How long would it take to reach that steady-state with time-stepping?\n", "\n", "We chan check by tracking the norm of `F(x,p)`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "Δt = 40000u\"yr\" |> u\"s\" |> ustrip\n", "x_hist, t_hist = time_steps(x₀, Δt, 4000, F, ∇ₓF)\n", "PyPlot.figure(figsize=(7,3))\n", "PyPlot.semilogy(t_hist .* 1u\"s\" .|> u\"yr\" .|> ustrip, [norm(F(s,p)) for i in 1:size(x_hist,2)], label=\"steady-state\")\n", "PyPlot.semilogy(t_hist .* 1u\"s\" .|> u\"yr\" .|> ustrip, [norm(F(x_hist[:,i],p)) for i in 1:size(x_hist,2)], label=\"time-stepping\")\n", "PyPlot.xlabel(\"simulation time (years)\"); PyPlot.ylabel(\"|F(x,p)| (mol m⁻³ s⁻¹)\");\n", "PyPlot.legend(); PyPlot.title(\"Stability of ¹⁴C age with simulation time\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Circulation 2: OCIM1\n", "\n", "The Ocean Circulation Inverse Model (OCIM) version 1 is loaded via" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "wet3D, grd, T = OCIM1.load() ;" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Redefine `F` and `∇ₓF` for the new `T`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "F, ∇ₓF = state_function_and_Jacobian(p -> T, G) ;" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Redefine `iwet` and `x` for the new grid size" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "iwet = findall(wet3D)\n", "x = p.Ratm * ones(length(iwet))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Redefine `h`, `z₀`, and `z` for the new grid" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "p.h = grd.δdepth[1] |> upreferred |> ustrip\n", "z = grd.depth_3D[iwet]\n", "z₀ = z[1]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Solve for the steady-state of radio carbon and convert to age" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "prob = SteadyStateProblem(F, ∇ₓF, x, p)\n", "s = solve(prob, CTKAlg()).u \n", "C14age = log.(p.Ratm ./ s) * p.τ * u\"s\" .|> u\"yr\"" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "And plot horizontal slices using GR and Interact:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "lon, lat = grd.lon |> ustrip, grd.lat |> ustrip\n", "function horizontal_slice(x, levels, title)\n", " x_3D = fill(NaN, size(grd))\n", " x_3D[iwet] .= x .|> ustrip\n", " GR.figure(figsize=(10,5))\n", " @manipulate for iz in 1:size(grd)[3]\n", " GR.xlim([0,360])\n", " GR.title(string(title, \" at $(AIBECS.round(grd.depth[iz])) depth\"))\n", " GR.contourf(lon, lat, x_3D[:,:,iz]', levels=levels)\n", " end\n", "end \n", "horizontal_slice(C14age, 0:100:2400, \"14C age [yr] using the OCIM1 circulation\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Or zonal slices:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "function zonal_slice(x, levels, title)\n", " x_3D = fill(NaN, size(grd))\n", " x_3D[iwet] .= x .|> ustrip\n", " GR.figure(figsize=(10,5))\n", " @manipulate for longitude in (grd.lon .|> ustrip)\n", " GR.title(string(title, \" at $(round(longitude))°\"))\n", " ilon = findfirst(ustrip.(grd.lon) .== longitude)\n", " GR.contourf(lat, reverse(-grd.depth .|> ustrip), reverse(x_3D[:,ilon,:], dims=2), levels=levels)\n", " end\n", "end \n", "zonal_slice(C14age, 0:100:2400, \"14C age [yr] using the OCIM1 circulation\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "

Example 2: A phosphorus cycle

\n", "\n", "Dissolved inorganic phosphrous (**DIP**)
\n", "(transported by the **ocean circulation**)\n", "\n", "$$\\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})},$$\n", "\n", "and particulate organic phosphrous (**POP**)
\n", "(transported by **sinking particles**)\n", "\n", "\n", "$$\\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})}.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Ocean circulation\n", "\n", "For DIP, the **advective–eddy-diffusive transport** operator, $\\nabla \\cdot (\\boldsymbol{u} + \\mathbf{K}\\cdot\\nabla)$, is converted into the matrix `T`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "T_DIP(p) = T" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Sinking particles\n", "\n", "For POP, the **particle flux divergence (PFD)** operator, $\\nabla \\cdot \\boldsymbol{w}$, is created via `buildPFD`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "T_POP(p) = buildPFD(grd, wet3D, sinking_speed = w(p))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The settling velocity, `w(p)`, is assumed linearly increasing with depth `z` to yield a \"Martin curve profile\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "w(p) = p.w₀ .+ p.w′ * (z .|> ustrip)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Uptake:\n", "$$\\frac{1}{τ} \\, \\frac{\\boldsymbol{x}_\\mathrm{DIP}^2}{\\boldsymbol{x}_\\mathrm{DIP} + k}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "relu(x) = (x .≥ 0) .* x\n", "zₑ = 80u\"m\" # depth of the euphotic zone\n", "function uptake(DIP, p)\n", " τ, k = p.τ, p.k\n", " DIP⁺ = relu(DIP)\n", " return 1/τ * DIP⁺.^2 ./ (DIP⁺ .+ k) .* (z .≤ zₑ)\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Remineralization:\n", "$$\\kappa \\, \\boldsymbol{x}_\\mathrm{POP}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "remineralization(POP, p) = p.κ * POP" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### \"Geological\" restoring\n", "$$\\frac{x_\\mathrm{geo} - \\boldsymbol{x}_\\mathrm{DIP}}{\\tau_\\mathrm{geo}}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "geores(x, p) = (p.xgeo .- x) / p.τgeo" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Net sum of local sources and sinks" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "G_DIP(DIP, POP, p) = -uptake(DIP, p) + remineralization(POP, p) + geores(DIP, p)\n", "G_POP(DIP, POP, p) = uptake(DIP, p) - remineralization(POP, p)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "t = empty_parameter_table() # empty table of parameters\n", "add_parameter!(t, :xgeo, 2.12u\"mmol/m^3\")\n", "add_parameter!(t, :τgeo, 1.0u\"Myr\")\n", "add_parameter!(t, :k, 6.62u\"μmol/m^3\")\n", "add_parameter!(t, :w₀, 0.64u\"m/d\")\n", "add_parameter!(t, :w′, 0.13u\"1/d\")\n", "add_parameter!(t, :κ, 0.19u\"1/d\")\n", "add_parameter!(t, :τ, 236.52u\"d\")\n", "initialize_Parameters_type(t, \"Pcycle_Parameters\") # Generate the parameter type\n", "p = Pcycle_Parameters()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# State function `F` and Jacobian `∇ₓF`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "nb = length(iwet); x = ones(2nb)\n", "F, ∇ₓF = state_function_and_Jacobian((T_DIP,T_POP), (G_DIP,G_POP), nb) ;" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Solve the steady-state PDE system" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "prob = SteadyStateProblem(F, ∇ₓF, x, p)\n", "s = solve(prob, CTKAlg(), preprint=\" \").u" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Plot DIP" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "DIP, POP = state_to_tracers(s, nb, 2)\n", "DIP_unit = u\"mmol/m^3\"\n", "horizontal_slice(DIP * u\"mol/m^3\" .|> DIP_unit, 0:0.3:3, \"P-cycle model: DIP [mmol m^{-3}]\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Plot POP" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "zonal_slice(POP .|> relu .|> log10, -10:1:10, \"P-cycle model: POP [log\\\\(mol m^{-3}\\\\)]\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We can also make publication-quality plots using Python's Cartopy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "ccrs = pyimport(\"cartopy.crs\")\n", "cfeature = pyimport(\"cartopy.feature\")\n", "lon_cyc = [lon; 360+lon[1]]\n", "DIP_3D = fill(NaN, size(grd)); DIP_3D[iwet] = DIP * u\"mol/m^3\" .|> u\"mmol/m^3\" .|> ustrip\n", "DIP_3D_cyc = cat(DIP_3D, DIP_3D[:,[1],:], dims=2)\n", "f1 = PyPlot.figure(figsize=(12,5))\n", "@manipulate for iz in 1:24\n", " withfig(f1, clear=true) do\n", " ax = PyPlot.subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))\n", " ax.add_feature(cfeature.COASTLINE, edgecolor=\"#000000\") # black coast lines\n", " ax.add_feature(cfeature.LAND, facecolor=\"#CCCCCC\") # gray land\n", " plt = PyPlot.contourf(lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=0:0.1:3.5, \n", " transform=ccrs.PlateCarree(), zorder=-1, extend=\"both\")\n", " plt2 = PyPlot.contour(plt, lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=plt.levels[1:5:end], \n", " transform=ccrs.PlateCarree(), zorder=0, extend=\"both\", colors=\"k\")\n", " ax.clabel(plt2, fmt=\"%2.1f\", colors=\"k\", fontsize=14)\n", " cbar = PyPlot.colorbar(plt, orientation=\"vertical\", extend=\"both\", ticks=plt2.levels)\n", " cbar.add_lines(plt2)\n", " cbar.ax.tick_params(labelsize=14)\n", " cbar.set_label(label=\"mmol / m³\", fontsize=16)\n", " PyPlot.title(\"Publication-quality DIP with Cartopy; depth = $(string(round(typeof(1u\"m\"),grd.depth[iz])))\", fontsize=20)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.1.1", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.1" }, "rise": { "scroll": true } }, "nbformat": 4, "nbformat_minor": 4 }