River discharge

In this tutorial we will simulate a fictitious radioactive tracer that is injected into the ocean by the 200 largest rivers (by estimated discharge). The 200 major rivers dataset from Dai and Trenberth (2002) is available from within the AIBECS. Once "born", our ficitious tracer decays with a parameter timescale $\tau$ as it flows through ocean basins.

The 3D tracer equation is:

\[\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \mathbf{K}\nabla)\right] x = s_\mathsf{rivers} - x / \tau\]

where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ represents the ocean circulation transport. (Tracer transport operators are described in the documentation.) The riverine source of the tracer is $s_\mathsf{rivers}$, and $x / \tau$ is the decay rate.

In AIBECS, we must recast this equation in the generic form

\[\left[\frac{\partial}{\partial t} + \mathbf{T}(\boldsymbol{p})\right] \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x},\boldsymbol{p}).\]

We start by telling Julia that we want to use the AIBECS and the OCIM2 transport matrix for the ocean circulation.

using AIBECS
grd, T_OCIM2 = OCIM2.load()

(, sparse([1, 2, 10384, 10442, 10443, 20825, 20883, 1, 2, 3  …  200160, 197886, 199766, 199777, 199778, 199779, 199790, 200156, 200159, 200160], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  200159, 200160, 200160, 200160, 200160, 200160, 200160, 200160, 200160, 200160], [0.00019778421518954799, 2.3427916742722093e-9, -1.9599474163829085e-7, -0.00019161212648881556, 4.8096149072091506e-9, -1.830592653460076e-9, 5.007679174162751e-9, -5.025164843241415e-8, 0.00018753126417941492, 4.264266869682882e-8  …  -2.196560075226544e-8, 1.0819937104262028e-10, 6.709812718407374e-9, -1.263521554746615e-9, -3.3927920410468295e-9, 7.593163378667893e-9, -7.410175543096161e-9, -3.441057669604186e-8, -2.0030251520181335e-8, 5.2794476107904204e-8], 200160, 200160))

For the radioactive decay, we simply use

function decay(x, p)
    @unpack τ = p
    return x / τ
end

decay (generic function with 1 method)

To build the river sources, we will load the geographic locations and discharge (in m³ s⁻¹) from the Dai and Trenberth (2017) dataset.

RIVERS = Rivers.load()

200-element Vector{AIBECS.Rivers.River{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, s⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}}:
 Amazon (-2.0, -55.5) 210472 m³ s⁻¹
 Congo (-4.3, 15.3) 41448 m³ s⁻¹
 Orinoco (8.1, -63.6) 35776 m³ s⁻¹
 Changjiang (30.8, 117.6) 29914 m³ s⁻¹
 Brahmaputra (25.2, 89.7) 19900 m³ s⁻¹
 Mississippi (32.3, -90.9) 19330 m³ s⁻¹
 Yenisey (67.4, 86.5) 18981 m³ s⁻¹
 Paraná (-32.7, -60.7) 17999 m³ s⁻¹
 Lena (70.7, 127.4) 16826 m³ s⁻¹
 Mekong (15.1, 105.8) 16636 m³ s⁻¹
 ⋮
 Guadalquivir (37.5, -6.0) 105 m³ s⁻¹
 Yaqui (29.2, -109.5) 105 m³ s⁻¹
 Colorado-TX (29.3, -96.1) 82 m³ s⁻¹
 Rio Grande (25.9, -97.4) 48 m³ s⁻¹
 de Grey (-20.3, 119.2) 44 m³ s⁻¹
 Gascoyne (-24.8, 113.8) 22 m³ s⁻¹
 Avon (-31.8, 116.1) 22 m³ s⁻¹
 Fortescue (-21.3, 116.2) 10 m³ s⁻¹
 Murchison (-27.9, 114.5) 6 m³ s⁻¹

This is an array of rivers, for which the type River{T} contains the river's name, lat–lon coordinates, and discharge in m³ s⁻¹. For example, the first river is the Amazon

r = RIVERS[1]

Amazon (-2.0, -55.5) 210472 m³ s⁻¹

We can check the locations with

using Plots
scatter([r.lon for r in RIVERS], [r.lat for r in RIVERS],
        zcolor=log10.(uconvert.(NoUnits, [r.VFR for r in RIVERS] / u"m^3/s")),
        colorbartitle="log₁₀(discharge / (1 m³ s⁻¹))")

We can regrid these into the OCIM2 grid and return the corresponding vector with

rivers = regrid(RIVERS, grd)

200160-element Vector{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, s⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}:
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
          ⋮
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹
 0.0 m³ s⁻¹

(Note this regridding uses NearestNeighbors.jl to assign a wet box as the mouth of each river, which sometimes is not exactly the real location of the river mouth.)

We control the global magnitude of the river discharge, $\sigma$ (in mol s⁻¹), by making it a parameter of our model. For that, we separate the river source

\[s_\mathsf{rivers} = \sigma s_0\]

into global magnitude ($\sigma$) and spatial pattern ($s_0$).

Since $\int s_0 \mathrm{d}V = 1$, s_0 can be computed by normalizing rivers. In Julia/AIBECS, this can be done by dividing rivers by the dot product v ⋅ rivers (or v'rivers in matrix form). (v ⋅ x is the discrete equivalent of the volume integral $\int x \mathrm{d}V$.)

v = vector_of_volumes(grd)
s_0 = rivers / (v'rivers)
function s_rivers(p)
    @unpack σ = p
    return σ * ustrip.(s_0) # we must remove the units in AIBECS here :(
end

s_rivers (generic function with 1 method)

We then write the generic $\boldsymbol{G}$ function, which is

G_radiorivers(x,p) = s_rivers(p) - decay(x,p)

G_radiorivers (generic function with 1 method)

Parameters

We specify some initial values for the parameters and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct RadioRiversParameters{U} <: AbstractParameters{U}
    τ::U | 5.0 | u"yr"
    σ::U | 1.0 | u"Gmol/yr"
end

initial_value (generic function with 51 methods)

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

p = RadioRiversParameters()

 Row │ Symbol  Value    Initial value  Unit
     │ Symbol  Float64  Float64        FreeUnit…
─────┼───────────────────────────────────────────
   1 │ τ           5.0            5.0  yr
   2 │ σ           1.0            1.0  Gmol yr⁻¹

We generate the state function F that gives the tendencies,

F = AIBECSFunction(T_OCIM2, G_radiorivers)

(::SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, AIBECS.var"#f#58"{Tuple{AIBECS.var"#49#50"{SparseMatrixCSC{Float64, Int64}}}, Vector{Int64}, AIBECS.var"#G#56"{Tuple{typeof(Main.G_radiorivers)}, 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{AIBECS.var"#49#50"{SparseMatrixCSC{Float64, Int64}}}, Int64, Vector{Int64}}, AIBECS.var"#∇ₓG#59"{Tuple{typeof(Main.G_radiorivers)}, 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 prob,

nb = sum(iswet(grd))
x = ones(nb) # initial guess
prob = SteadyStateProblem(F, x, p)

SteadyStateProblem with uType Vector{Float64}. In-place: false
u0: 200160-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 ⋮
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

and solve it

sol = solve(prob, CTKAlg()).u * u"mol/m^3"

200160-element Vector{Quantity{Float64, 𝐍 𝐋⁻³, Unitful.FreeUnits{(m⁻³, mol), 𝐍 𝐋⁻³, nothing}}}:
  -3.597832067855361e-12 mol m⁻³
  -5.353871365614832e-12 mol m⁻³
 -6.4834851941454704e-12 mol m⁻³
  -7.823908119438312e-12 mol m⁻³
   -9.79044237114253e-12 mol m⁻³
  -9.770934123554808e-12 mol m⁻³
 -2.4448173185331663e-11 mol m⁻³
  -6.738392003832571e-11 mol m⁻³
   -6.51083745081848e-11 mol m⁻³
  -7.890402312090086e-11 mol m⁻³
                               ⋮
  -3.876074162661911e-13 mol m⁻³
 -1.3409356343223646e-13 mol m⁻³
  -7.257906578905357e-13 mol m⁻³
  -6.182255119626855e-13 mol m⁻³
 -3.4001872611692876e-13 mol m⁻³
    6.19239679441887e-13 mol m⁻³
  -5.257895087043751e-13 mol m⁻³
  -6.838780987222902e-13 mol m⁻³
  -4.278816669049942e-13 mol m⁻³

Let's now run some visualizations using the plot recipes. Taking a horizontal slice of the 3D field at 200m gives

cmap = :viridis
plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=200, color=cmap)

and at 500m:

plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,0.05))

Or we can change the timescale and watch the tracer fill the oceans:

p = RadioRiversParameters(τ = 50.0u"yr")
prob = SteadyStateProblem(F, x, p)
sol_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,1))

Point-wise sources like the rivers used here can be problematic numerically because these can generate large gradients and numerical noise with the given spatial discretization:

cmap = :RdYlBu_4
plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))

Note, for an example of numerical noise, the negative values appearing in red. Hence, it might be wise to smooth the 3D field of the river sources. We can do this using the smooth_operator function (smooth_operator creates an operator that can be applied to smooth the AIBECS vector along the stencil of the given transport matrix and in a volume-conservative way). And we can check that the negative values disappear after smoothing of the sources:

S = smooth_operator(grd, T_OCIM2)
s_0 = S * (S * s_0) # smooth river sources twice
sol_smooth = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))

Note that such numerical noise generally dampens out with distance and depth

cmap = :viridis
plot(
    plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
    plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
    layout=(2,1),
    size=(600,600)
)

Nevertheless, it can be good to reduce noise as much as possible!

This page was generated using Literate.jl.