A dust model

In this tutorial we will simulate particles of dust that graviationally settle. 2D maps of dust deposition are available from within the AIBECS. To use these, we will convert them into sources in the top layer of the model grid. Once "born", depending on the settling velocity, these dust particles can spread horizontally until they eventually reach the sea floor and are buried in the sediments.

The governing equation of the 3D concentration field of dust, denoted , is:

where and represent the ocean circulation transport and the vertical settling of dust particles, respectively. (Tracer transport operators are described in the documentation.) The only source term, , is the dust deposition.

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

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

julia

using AIBECS
using JLD2       # required by `OCIM0.load`
using NCDatasets # required by `AeolianSources.load`
grd, T_OCIM = OCIM0.load()

(, sparse([1, 2, 57, 10172, 10229, 1, 2, 3, 58, 10173  …  191168, 190655, 191165, 191167, 191168, 191169, 190656, 191166, 191168, 191169], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2  …  191167, 191168, 191168, 191168, 191168, 191168, 191169, 191169, 191169, 191169], [3.4355389417142345e-7, 6.157619696306188e-8, -2.61541451353055e-7, -7.123747549591164e-7, 5.461579111836965e-7, -1.0940833534352636e-7, 3.1285031065724247e-7, 1.9807440791995236e-9, -2.4207715692543693e-8, -2.3441757104252792e-7  …  -1.815611679276709e-7, -2.3262536747686192e-8, -1.0813061941848588e-7, 1.421122420594373e-7, 6.118547004578337e-8, -7.254164692789022e-8, -1.8087448238776243e-9, -7.125518535170859e-8, 3.2353843704816876e-8, 4.072334915760755e-8], 191169, 191169))

The transport of dust is the sum of the ocean circulation and a settling-transport:

julia

function T_dust(p)
    @unpack fsedremin = p
    return LinearOperators((T_OCIM, transportoperator(grd, z -> w(z, p), fsedremin = fsedremin)))
end

T_dust (generic function with 1 method)

For the settling transport of dust particles, we have used the transportoperator function for which we need to define the settling velocity w(z,p) as a function of depth z and of the parameters p.

Note

We have imposed a fsedremin remineralization at the sediment with its keyword argument. (Default is fsedremin=1.0.) This turns the settling transport into a sink that removes all the dust reaching the sea floor.

Following the assumption that increases linearly with depth, we write it as

julia

function w(z, p)
    @unpack w₀, w′ = p
    return w₀ + w′ * z
end

w (generic function with 1 method)

The parameters fsedremin, w₀, and w′ will be defined shortly.

Dust depositon

The AIBECS allows you to use climatological dust deposition fields (from the Chien et al. dataset by default).

julia

s_A_2D = AeolianSources.load()

Dict{Symbol, Array{Float64}} with 9 entries:
  :dust => [3.47645e-14 4.67411e-14 … 2.70313e-12 1.92588e-12; 3.47699e-14 4.56…
  :fire => [5.08321e-16 7.00633e-16 … 6.73572e-14 6.07252e-14; 5.08301e-16 6.67…
  :volc => [6.72953e-16 9.75435e-16 … 1.32062e-14 8.07187e-15; 6.72933e-16 8.97…
  :lat  => [-90.0, -88.1053, -86.2105, -84.3158, -82.4211, -80.5263, -78.6316, …
  :bioF => [6.21323e-17 1.00159e-16 … 5.83582e-15 4.04055e-15; 6.2127e-17 9.280…
  :salt => [8.30242e-13 1.30496e-12 … 3.43202e-11 2.11114e-11; 8.30477e-13 1.22…
  :FF   => [2.01311e-16 2.93369e-16 … 4.19621e-14 2.78746e-14; 2.01282e-16 2.73…
  :lon  => [0.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20.0, 22.5  …  335.0, 3…
  :bioP => [7.59786e-16 1.0005e-15 … 1.75302e-13 1.51327e-13; 7.59332e-16 9.406…

s_A_2D[:dust] is a 2D annual-mean map (kg m⁻² s⁻¹). We transpose it to (lat, lon) order and regrid it onto the surface layer of grd, the 3D grid of the ocean circulation:

julia

s_dust_2D = regrid(permutedims(s_A_2D[:dust], (2, 1)), s_A_2D[:lat], s_A_2D[:lon], grd)

90×180 Matrix{Float64}:
 4.0861e-14   4.02723e-14  3.92638e-14  …  4.23192e-14  4.14034e-14
 5.58523e-14  5.28017e-14  4.9783e-14      6.29739e-14  5.93423e-14
 5.58088e-14  5.2256e-14   4.86971e-14     6.61204e-14  6.0421e-14
 4.4301e-14   3.99454e-14  3.605e-14       4.99711e-14  4.73808e-14
 3.14211e-14  2.86422e-14  2.65818e-14     3.76665e-14  3.42194e-14
 3.3941e-14   3.02657e-14  2.6606e-14   …  4.51644e-14  3.8863e-14
 4.53696e-14  4.2776e-14   3.72113e-14     5.3644e-14   4.83879e-14
 5.31449e-14  4.80786e-14  4.30237e-14     5.9285e-14   5.59273e-14
 1.4441e-13   1.42862e-13  1.41995e-13     1.67321e-13  1.51856e-13
 1.34627e-13  1.33476e-13  1.34818e-13     1.66875e-13  1.4447e-13
 ⋮                                      ⋱               
 3.74026e-12  4.02252e-12  4.5013e-12      3.60454e-12  3.61945e-12
 4.24891e-12  4.45213e-12  4.66692e-12     3.80302e-12  4.02533e-12
 4.11568e-12  4.54572e-12  4.8539e-12      3.36415e-12  3.72062e-12
 3.47995e-12  3.84901e-12  3.92341e-12     2.71971e-12  3.05456e-12
 2.67713e-12  2.84816e-12  2.96605e-12  …  2.28933e-12  2.47973e-12
 2.26981e-12  2.34474e-12  2.41625e-12     2.07681e-12  2.17852e-12
 2.17696e-12  2.22869e-12  2.27043e-12     2.09304e-12  2.13253e-12
 2.46303e-12  2.4854e-12   2.46811e-12     2.3863e-12   2.42223e-12
 2.3294e-12   2.32015e-12  2.32332e-12     2.29825e-12  2.32053e-12

Then we create a blank 3D array of zeros, of which we paint the top layer:

julia

s_dust_3D = zeros(size(grd)...)
s_dust_3D[:, :, 1] .= s_dust_2D

90×180 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 4.0861e-14   4.02723e-14  3.92638e-14  …  4.23192e-14  4.14034e-14
 5.58523e-14  5.28017e-14  4.9783e-14      6.29739e-14  5.93423e-14
 5.58088e-14  5.2256e-14   4.86971e-14     6.61204e-14  6.0421e-14
 4.4301e-14   3.99454e-14  3.605e-14       4.99711e-14  4.73808e-14
 3.14211e-14  2.86422e-14  2.65818e-14     3.76665e-14  3.42194e-14
 3.3941e-14   3.02657e-14  2.6606e-14   …  4.51644e-14  3.8863e-14
 4.53696e-14  4.2776e-14   3.72113e-14     5.3644e-14   4.83879e-14
 5.31449e-14  4.80786e-14  4.30237e-14     5.9285e-14   5.59273e-14
 1.4441e-13   1.42862e-13  1.41995e-13     1.67321e-13  1.51856e-13
 1.34627e-13  1.33476e-13  1.34818e-13     1.66875e-13  1.4447e-13
 ⋮                                      ⋱               
 3.74026e-12  4.02252e-12  4.5013e-12      3.60454e-12  3.61945e-12
 4.24891e-12  4.45213e-12  4.66692e-12     3.80302e-12  4.02533e-12
 4.11568e-12  4.54572e-12  4.8539e-12      3.36415e-12  3.72062e-12
 3.47995e-12  3.84901e-12  3.92341e-12     2.71971e-12  3.05456e-12
 2.67713e-12  2.84816e-12  2.96605e-12  …  2.28933e-12  2.47973e-12
 2.26981e-12  2.34474e-12  2.41625e-12     2.07681e-12  2.17852e-12
 2.17696e-12  2.22869e-12  2.27043e-12     2.09304e-12  2.13253e-12
 2.46303e-12  2.4854e-12   2.46811e-12     2.3863e-12   2.42223e-12
 2.3294e-12   2.32015e-12  2.32332e-12     2.29825e-12  2.32053e-12

we convert it to a SI volumetric unit (i.e., in g m⁻³)

julia

s_dust_3D = ustrip.(upreferred.(s_dust_3D * u"kg/m^2/s" ./ grd.δz_3D))

90×180×24 Array{Float64, 3}:
[:, :, 1] =
 1.13078e-15  1.11449e-15  1.08658e-15  …  1.17114e-15  1.14579e-15
 1.54565e-15  1.46123e-15  1.37769e-15     1.74273e-15  1.64223e-15
 1.54445e-15  1.44613e-15  1.34764e-15     1.82981e-15  1.67208e-15
 1.22598e-15  1.10544e-15  9.97645e-16     1.3829e-15   1.31121e-15
 8.69544e-16  7.92642e-16  7.35623e-16     1.04238e-15  9.46985e-16
 9.3928e-16   8.37571e-16  7.36293e-16  …  1.24988e-15  1.07549e-15
 1.25556e-15  1.18378e-15  1.02978e-15     1.48454e-15  1.33908e-15
 1.47073e-15  1.33052e-15  1.19063e-15     1.64065e-15  1.54773e-15
 3.99639e-15  3.95355e-15  3.92957e-15     4.63042e-15  4.20245e-15
 3.72566e-15  3.69379e-15  3.73094e-15     4.61808e-15  3.99804e-15
 ⋮                                      ⋱               
 1.03508e-13  1.11319e-13  1.24568e-13     9.97516e-14  1.00164e-13
 1.17584e-13  1.23208e-13  1.29152e-13     1.05244e-13  1.11397e-13
 1.13897e-13  1.25798e-13  1.34327e-13     9.30992e-14  1.02964e-13
 9.63039e-14  1.06517e-13  1.08576e-13     7.5265e-14   8.45318e-14
 7.40866e-14  7.88196e-14  8.20823e-14  …  6.33547e-14  6.86238e-14
 6.28146e-14  6.48882e-14  6.68671e-14     5.74734e-14  6.02881e-14
 6.02451e-14  6.16767e-14  6.28316e-14     5.79227e-14  5.90154e-14
 6.81616e-14  6.87808e-14  6.83021e-14     6.60383e-14  6.70327e-14
 6.44636e-14  6.42077e-14  6.42953e-14     6.36017e-14  6.42183e-14

[:, :, 2] =
 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  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     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  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  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  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  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  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  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  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  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  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  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  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  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

[:, :, 3] =
 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  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     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  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  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  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  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  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  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  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  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  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  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  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  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

;;; … 

[:, :, 22] =
 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  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     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  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  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  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  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  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  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  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  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  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  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  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  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

[:, :, 23] =
 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  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     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  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  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  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  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  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  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  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  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  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  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  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  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

[:, :, 24] =
 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  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     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  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  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  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  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  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  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  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  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  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  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  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  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

Finally, we vectorize it via

julia

s_dust = vectorize(s_dust_3D, grd) # which is the same as `s_dust_3D[iswet(grd)]`

191169-element Vector{Float64}:
 6.515217379785785e-15
 1.0037390656128683e-14
 1.3316895771101413e-14
 1.804431065139669e-14
 2.370608647844632e-14
 3.5478145460585265e-14
 4.647631247106428e-14
 4.5066419096689556e-14
 4.729787372474041e-14
 4.912364229160753e-14
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

We then write the generic function, which is simply the constant dust source:

julia

G_dust(x, p) = s_dust

G_dust (generic function with 1 method)

Parameters

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

julia

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct DustModelParameters{U} <: AbstractParameters{U}
    w₀::U | 0.1 | u"km/yr"
    w′::U | 0.1 | u"km/yr/km"
    fsedremin::U | 5.0 | u"percent"
end

initial_value (generic function with 27 methods)

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

julia

p = DustModelParameters()

Main.DustModelParameters{Float64}
  w₀ = 0.1 (km yr⁻¹)
  w′ = 0.1 (yr⁻¹)
  fsedremin = 5.0 (%)

We build the state function F,

julia

nb = count(iswet(grd))
F = AIBECSFunction(T_dust, G_dust, nb)

SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, AIBECS.var"#f#71"{Tuple{typeof(Main.T_dust)}, Vector{Int64}, AIBECS.var"#G#68"{Tuple{typeof(Main.G_dust)}, AIBECS.var"#tracers#66"{Int64, Int64}}, AIBECS.var"#tracer#67"{Int64, Int64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, AIBECS.var"#jac#81"{AIBECS.var"#T#76"{Tuple{typeof(Main.T_dust)}, Int64, Vector{Int64}}, AIBECS.var"#72#73"{Tuple{typeof(Main.G_dust)}, Int64, Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}(AIBECS.var"#f#71"{Tuple{typeof(Main.T_dust)}, Vector{Int64}, AIBECS.var"#G#68"{Tuple{typeof(Main.G_dust)}, AIBECS.var"#tracers#66"{Int64, Int64}}, AIBECS.var"#tracer#67"{Int64, Int64}}((Main.T_dust,), [1], AIBECS.var"#G#68"{Tuple{typeof(Main.G_dust)}, AIBECS.var"#tracers#66"{Int64, Int64}}((Main.G_dust,), AIBECS.var"#tracers#66"{Int64, Int64}(191169, 1)), AIBECS.var"#tracer#67"{Int64, Int64}(191169, 1)), LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, AIBECS.var"#jac#81"{AIBECS.var"#T#76"{Tuple{typeof(Main.T_dust)}, Int64, Vector{Int64}}, AIBECS.var"#72#73"{Tuple{typeof(Main.G_dust)}, Int64, Int64}}(AIBECS.var"#T#76"{Tuple{typeof(Main.T_dust)}, Int64, Vector{Int64}}((Main.T_dust,), 1, [1]), AIBECS.var"#72#73"{Tuple{typeof(Main.G_dust)}, Int64, Int64}((Main.G_dust,), 191169, 1)), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing, nothing, nothing)

the steady-state problem,

julia

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

SteadyStateProblem with uType Vector{Float64}. In-place: false
u0: 191169-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

julia

sol = solve(prob, CTKAlg()).u

191169-element Vector{Float64}:
 1.2200651489102654e-7
 1.5330849345573938e-7
 1.7141350938665312e-7
 1.5349199999903362e-7
 2.725786698576371e-7
 3.856189246809194e-7
 4.90324899460326e-7
 5.89799068806407e-7
 6.009599564649193e-7
 7.529330591227179e-7
 ⋮
 1.1485079759454564e-5
 1.4445383317851501e-5
 1.46217710297325e-5
 6.9338103768634696e-6
 7.1223364274120016e-6
 7.008123986240733e-6
 7.063844599031199e-6
 6.9683554916322166e-6
 7.038666482883649e-6

Let's now run some vizualizations using the Plots.jl recipes.

julia

using Plots

Let's have a look at a map of the mean dust concentration:

julia

plotverticalmean(sol * u"g/m^3", grd, color = cgrad(:turbid, rev = true, scale = :exp))

compared to the original source

julia

plotverticalintegral(s_dust * u"g/m^3/s", xlabel = "", ylabel = "", grd, color = cgrad(:starrynight, scale = :exp), title = "Surface flux", zunit = u"mg/m^2/yr")

Huh? This is odd isn't it? What is happening? It seems that dust settles slowly enough to actually explore the ocean basins a little bit. This is of course dependent on our choice of model parameters.

Let's look at what is exported below 500 m to compare:

julia

wflux = map(z -> w(z, p), depthvec(grd)) * u"m/s"
plothorizontalslice(sol * u"g/m^3" .* wflux, grd, depth = 500u"m", color = cgrad(:starrynight, scale = :exp), xlabel = "", ylabel = "", title = "Flux at 500m", zunit = u"mg/yr/m^2")

We can also take a look at the global mean profile (the horizontal average)

julia

profile_plot = plothorizontalmean(sol * u"g/m^3" .|> u"mg/m^3", grd)

As you may see, it also seems that with this model, dust (here with a 1% fraction allowed to stay in the system when it reaches the seafloor) accumulates a little bit in the deepest boxes of the model grid.

Don't hesitate to play around with the parameters and see how things change! For example, impose a larger fraction that stays in the bottom or change other parameters via

julia

p2 = DustModelParameters(w₀ = 0.1u"km/yr", w′ = 0.1u"km/yr/km", fsedremin = 95.0u"percent")
prob2 = SteadyStateProblem(F, x, p2)
sol2 = solve(prob2, CTKAlg()).u
plotverticalmean(sol2 * u"g/m^3", grd, color = cgrad(:turbid, rev = true, scale = :exp))

And watch how all the dust accumulates in the deepest holes of the ocean basins!

Indeed, this dust tracer penetrates much deeper in the ocean:

julia

plothorizontalmean!(profile_plot, sol2 * u"g/m^3" .|> u"mg/m^3", grd)

This page was generated using Literate.jl.

A dust model ​

Dust depositon ​

Parameters ​

A dust model

Dust depositon

Parameters