Edit on GitHub

A dust model

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.)

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 $x_\mathsf{dust}$, is:

\[\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \boldsymbol{w} + \mathbf{K}\nabla)\right] x_\mathsf{dust} = s_\mathsf{dust}\]

where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ and $\nabla \cdot \boldsymbol{w}$ 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, $s_\mathsf{dust}$, is the dust deposition.

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 we want to use the AIBECS and the OCIM0.1 transport matrix for the ocean circulation.

using AIBECS
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:

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 $w(z) = w_0 + w' z$ increases linearly with depth, we write it as

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 Mahowald lab.

s_A_2D = AeolianSources.load()
Dict{Symbol, Array{Float64}} with 9 entries:
  :dust => [8.77097e-16 1.16976e-15 … 1.42237e-12 1.01189e-12; 8.71268e-16 1.12…
  :fire => [2.52328e-17 5.75708e-17 … 6.04506e-16 1.57569e-16; 2.51184e-17 5.38…
  :volc => [8.11113e-18 1.92487e-17 … 3.27646e-15 2.57447e-15; 8.06487e-18 1.79…
  :lat  => [-90.0, -88.1053, -86.2105, -84.3158, -82.4211, -80.5263, -78.6316, …
  :bioF => [6.68909e-18 1.52997e-17 … 3.27548e-15 2.10729e-15; 6.65695e-18 1.43…
  :salt => [4.40596e-13 5.79443e-13 … 1.83619e-11 9.72046e-12; 4.37885e-13 5.75…
  :FF   => [3.44167e-17 7.07536e-17 … 2.465e-14 1.64434e-14; 3.4266e-17 6.68044…
  :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 => [2.69616e-15 4.09648e-15 … 1.1505e-12 1.02498e-12; 2.67617e-15 3.944…

First, we must regrid these 2D maps onto the surface layer of grd, the 3D grid of the ocean circulation. To do that, we interpolate the 2D map of dust deposition onto the ocean grid's lat and lon:

s_dust_2D_monthly = s_A_2D[:dust] # kg m⁻² s⁻¹
s_dust_2D_annual = permutedims(dropdims(sum(s_dust_2D_monthly, dims=3), dims=3), (2,1)) / 12
s_dust_2D = regrid(s_dust_2D_annual, 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

(This required a bit of work to reshape s_dust_2D_annual before using regrid.)

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

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⁻³)

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.3757e-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.6938e-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.79226e-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

s_dust = vectorize(s_dust_3D, grd) # which is the same as `s_dust_3D[iswet(grd)]`
191169-element Vector{Float64}:
 6.515217452663623e-15
 1.0037390136081203e-14
 1.3316895913883588e-14
 1.8044310609691748e-14
 2.370608487073876e-14
 3.547814503526649e-14
 4.6476312590025904e-14
 4.506642014822272e-14
 4.729787498081552e-14
 4.9123644835055455e-14
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

We then write the generic $\boldsymbol{G}$ function, which is simply the constant dust source:

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.

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 49 methods)

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

p = DustModelParameters()
 Row │ Symbol     Value    Initial value  Unit
     │ Symbol     Float64  Float64        FreeUnit…
─────┼──────────────────────────────────────────────
   1 │ w₀             0.1            0.1  km yr⁻¹
   2 │ w′             0.1            0.1  yr⁻¹
   3 │ fsedremin      5.0            5.0  %

We build the state function F,

nb = count(iswet(grd))
F = AIBECSFunction(T_dust, G_dust, nb)
(::SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, AIBECS.var"#f#58"{Tuple{typeof(Main.T_dust)}, Vector{Int64}, AIBECS.var"#G#56"{Tuple{typeof(Main.G_dust)}, 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{typeof(Main.T_dust)}, Int64, Vector{Int64}}, AIBECS.var"#∇ₓG#59"{Tuple{typeof(Main.G_dust)}, 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,

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

sol = solve(prob, CTKAlg()).u
191169-element Vector{Float64}:
 1.2200651544363123e-7
 1.5330849234670112e-7
 1.7141351057258837e-7
 1.5349199665977283e-7
 2.725786649478154e-7
 3.8561893186674983e-7
 4.903249051708795e-7
 5.897990764299553e-7
 6.009599609290184e-7
 7.52933061274408e-7
 ⋮
 1.1485079770939733e-5
 1.4445383282368463e-5
 1.4621770991696269e-5
 6.933810375744731e-6
 7.1223364272353766e-6
 7.008123985363697e-6
 7.06384459855251e-6
 6.968355490653631e-6
 7.038666482342047e-6

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

using Plots

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

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

compared to the original source

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:

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)

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

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:

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

This page was generated using Literate.jl.