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 $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)))
endT_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.
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
endw (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_2D90×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-12we 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.0Finally, 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.0We then write the generic $\boldsymbol{G}$ function, which is simply the constant dust source:
G_dust(x, p) = s_dustG_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"
endinitial_value (generic function with 27 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, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, 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.0and solve it
sol = solve(prob, CTKAlg()).u191169-element Vector{Float64}:
1.2200651544363118e-7
1.5330849234670112e-7
1.7141351057258831e-7
1.5349199665977283e-7
2.7257866494781546e-7
3.856189318667499e-7
4.903249051708794e-7
5.897990764299552e-7
6.009599609290184e-7
7.529330612744083e-7
⋮
1.1485079770939731e-5
1.4445383282368463e-5
1.462177099169627e-5
6.933810375744733e-6
7.122336427235378e-6
7.0081239853637e-6
7.0638445985525105e-6
6.9683554906536335e-6
7.038666482342048e-6Let's now run some vizualizations using the Plots.jl recipes.
using PlotsLet'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.