Edit on GitHub

Groundwater discharge

Note

All the AIBECS tutorials and how-to guides are available as Jupyter notebooks. You can view them with nbvieweror execute them online with binder by clicking on the badges above! (Note that binder can be slow to launch and its memory caps can be a problem when running.)

Note

This tutorial is very similar to the river discharge notebook

In this tutorial we will simulate a fictitious radioactive tracer that is injected into the ocean by groundwater discharge. The global groundwater discharge dataset for 40,000 coastal watersheds from Luijendijk et al. (2020) 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{gw} - 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 source of the tracer is $s_\mathsf{gw}$, 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 groundwater sources, we will load the geographic locations and discharge of groundwater (in m³ yr⁻¹) from the Luijendijk et al. (2020)) dataset.

gws = GroundWaters.load()
40409-element Vector{AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}}:
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.1794887415, 24.2257701314, 47358.8148595 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.3417438136, 23.8746131396, 185142.093516 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.4161096656, 23.6050233764, 171000.541515 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.4266133609, 23.4030494462, 206889.096486 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.7820677063, 22.6137148344, 667545.345579 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.8962644211, 22.6987844211, 72625.201875 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.8614759879, 22.5822032466, 911480.517723 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(36.6601448123, 21.8794834664, 472586.358228 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(37.2143679963, 19.334562086, 1.33262583591e6 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(37.3177380952, 19.3185952381, 28384.2888 m³ yr⁻¹)
 ⋮
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(144.157841392, 72.2055303019, 11.7123541475 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(144.236000631, 72.2728170091, 41.4176247705 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(145.145353536, 71.8562156581, 73.1538069047 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(128.065238136, 73.2051409179, 105.676353367 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(143.155736084, 75.5195182604, 93.4943801401 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(143.243027112, 75.6684283401, 48.7752847061 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(142.916389905, 75.6972027214, 69.4594679685 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(128.501004892, 73.1684350529, 87.5895493016 m³ yr⁻¹)
 AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(141.832103547, 75.6537345857, 45.9949603817 m³ yr⁻¹)

This is an array of groundwater sources, for which the type GroundWaterSource{T} contains the lat–lon coordinates and discharge in m³ yr⁻¹. For example, the first element is

gws[1]
AIBECS.GroundWaters.GroundWaterSource{Quantity{Float64, 𝐋³ 𝐓⁻¹, Unitful.FreeUnits{(m³, yr⁻¹), 𝐋³ 𝐓⁻¹, nothing}}}(35.1794887415, 24.2257701314, 47358.8148595 m³ yr⁻¹)

We can check the locations with

using Plots
scatter([x.lon for x in gws], [x.lat for x in gws],
        zcolor=log10.(ustrip.([x.VFR for x in gws] / u"m^3/s")),
        colorbartitle="log₁₀(discharge / (1 m³ s⁻¹))",
        clim=(0,10), fmt=:png)

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

gws_grd = ustrip.(u"m^3/s", regrid(gws, grd))
200160-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

(We convert years to seconds for AIBECS)

We can control the global magnitude of the discharge by chosing a concentration of our tracer in groundwater, $C_\mathsf{gw}$, in mol m⁻³. For that, we need the volumes of the grid boxes

v = volumevec(grd)
200160-element Vector{Float64}:
 5.693492696598497e11
 6.267656666812283e11
 6.834351352130837e11
 7.392901411788552e11
 7.942641211123751e11
 8.482915614828765e11
 9.013080767686609e11
 9.532504861865775e11
 1.004056888985702e12
 1.053666738215696e12
 ⋮
 2.7437236229729617e13
 2.6923849140323855e13
 2.83653451428806e13
 2.7917925846214992e13
 2.7437236229729617e13
 2.6923849140323855e13
 2.83653451428806e13
 2.7917925846214992e13
 2.7437236229729617e13

because the source is given by the mass (discharge × concentration) divided by the volume of each box:

function s_gw(p)
    @unpack C_gw = p
    return @. gws_grd * C_gw / v
end
s_gw (generic function with 1 method)

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

G_gw(x,p) = s_gw(p) - decay(x,p)
G_gw (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 GroundWatersParameters{U} <: AbstractParameters{U}
    τ::U    | 20.0 | u"yr"
    C_gw::U |  1.0 | u"mol/m^3"
end
initial_value (generic function with 53 methods)

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

p = GroundWatersParameters()
 Row │ Symbol  Value    Initial value  Unit
     │ Symbol  Float64  Float64        FreeUnit…
─────┼───────────────────────────────────────────
   1 │ τ          20.0           20.0  yr
   2 │ C_gw        1.0            1.0  mol m⁻³

We build the state function F,

F = AIBECSFunction(T_OCIM2, G_gw)
(::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_gw)}, 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_gw)}, Int64, Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}) (generic function with 1 method)

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}}}:
 -1.1976569114245692e-7 mol m⁻³
 -1.4479734380088998e-7 mol m⁻³
 -1.4879453161967887e-7 mol m⁻³
 -1.5817358327495586e-7 mol m⁻³
 -1.9520226202176615e-7 mol m⁻³
 -2.0677685877025549e-7 mol m⁻³
 -2.2443661038761358e-7 mol m⁻³
 -4.1123842884196016e-7 mol m⁻³
  -5.599692326820508e-7 mol m⁻³
  -4.425251364360959e-7 mol m⁻³
                              ⋮
   5.356179456771745e-8 mol m⁻³
   9.890403277283755e-8 mol m⁻³
 -1.0615715734081943e-7 mol m⁻³
  -3.777197030624588e-8 mol m⁻³
   7.586448212525744e-8 mol m⁻³
  1.8682291679180856e-7 mol m⁻³
   -6.14559168215307e-8 mol m⁻³
   -6.28359735212723e-8 mol m⁻³
    6.06816725082375e-8 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, clim=(0,100))

and at 500m:

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

Or we can increase the decay timescale (×10) and decrease the groundwater concentration (÷10) to get a different (more well-mixed) tracer distribution:

p = GroundWatersParameters(τ = 200.0u"yr", C_gw = 0.1u"mol/m^3")
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,25))

This page was generated using Literate.jl.