Optimization of a simple phosphorus cycling model

Tracer equations

We consider a simple model for the cycling of phosphorus with 2 state variables consisting of phosphate (PO₄) AKA dissolved inorganic phosphorus (DIP) and particulate organic phosphorus (POP). The dissolved phases are transported by advection and diffusion whereas the particulate phase sinks rapidly down the water column without any appreciable transport by the circulation.

The governing equations that couple the 3D concentration fields of DIP and POP, denoted $x_\mathsf{DIP}$ and $x_\mathsf{POP}$, respectively, are:

and

The 3D vector $\boldsymbol{u}$ is the fluid velocity and $\mathbf{K}$ is the eddy-diffusion tensor. Thus, $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is a differential operator that represents the transport by the ocean circulation. The 3D vector $\boldsymbol{w}$ is the particulate sinking velocity, so that the operator $\nabla \cdot \boldsymbol{w}$ acting on $x_\mathsf{POP}$ is the divergence of the flux of POP. We will assume that it increases linearly with depth, i.e., that

where $\boldsymbol{w}_0$ and $\boldsymbol{w}'$ are simple scalars, independent of location $\boldsymbol{r}$. ($z$ represents the depth of location $\boldsymbol{r}=(x,y,z)$.)

The function $U$ represents the biological uptake of DIP by phytoplankton, which we model as

where the timescale, $\tau$, sets the maximum specific DIP uptake rate. The first term in the product scales linearly with DIP, which is a nice property to have in general. It also makes it easier to reason about the extreme cases. However, here the uptake is not linear: the second term in the product is a hyperbolic term (after Monod or Michaelis-Menten) and reduces the slope of the uptake function when DIP is low. When DIP is very low, the slope is close to 0. But when DIP is very high, the slope is close to $1/τ$. The last term is to be understood as a condition, i.e., that the uptake happens only in the euphotic layer, where there is sufficient light for photosynthesis to occur. The depth $z_0$ is the depth of the euphotic layer.

The function $R$ defines the remineralization rate of POP, which converts POP back into DIP. For the remineralization, we simply use a linear rate constant, i.e.,

Translation into AIBECS code

Transport matrices

We start by telling Julia we want to use the AIBECS

using AIBECS

We load the circulation and grid of the OCIM0.1 via

wet3D, grd, T_OCIM = OCIM0.load() ;

(Bool[false false … false false; false false … false false; … ; true true … true true; false false … false false]

Bool[false false … false false; false false … false false; … ; true true … true true; false false … false false]

Bool[false false … false false; false false … false false; … ; true true … true true; false false … false false]

...

Bool[false false … false false; false false … false false; … ; false false … false false; false false … false false]

Bool[false false … false false; false false … false false; … ; false false … false false; false false … false false]

Bool[false false … false false; false false … false false; … ; false false … false false; false false … false false], , 
  [1     ,      1]  =  3.43554e-7
  [2     ,      1]  =  6.15762e-8
  [57    ,      1]  =  -2.61541e-7
  [10172 ,      1]  =  -7.12375e-7
  [10229 ,      1]  =  5.46158e-7
  [1     ,      2]  =  -1.09408e-7
  [2     ,      2]  =  3.1285e-7
  [3     ,      2]  =  1.98074e-9
  [58    ,      2]  =  -2.42077e-8
  ⋮
  [190655, 191168]  =  -2.32625e-8
  [191165, 191168]  =  -1.08131e-7
  [191167, 191168]  =  1.42112e-7
  [191168, 191168]  =  6.11855e-8
  [191169, 191168]  =  -7.25416e-8
  [190656, 191169]  =  -1.80874e-9
  [191166, 191169]  =  -7.12552e-8
  [191168, 191169]  =  3.23538e-8
  [191169, 191169]  =  4.07233e-8)

where the operator $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is represented by the sparse matrix T_OCIM. In AIBECS, all the transport matrices must be declared as functions of the parameters p, so we simply define

T_DIP(p) = T_OCIM

T_DIP (generic function with 1 method)

For the sinking of particles, we are going to create a sparse matrix that depends on the parameters, too, which define how fast particles sink.

T_POP(p) = buildPFD(grd, wet3D, sinking_speed = w(p))

T_POP (generic function with 1 method)

for which we need to define the sinking speed w(p) as a function of the parameters p. Following the assumption that it increases linearly with depth, we write it as

w(p) = p.w₀ .+ p.w′ * z

w (generic function with 1 method)

For this to work, we must create a vector of depths, z, which is simply done via

iwet = findall(wet3D)
z = ustrip.(grd.depth_3D[iwet])

191169-element Array{Float64,1}:
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
   18.0675569520817
    ⋮              
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175
 5433.2531421838175

(We must strip the units via ustrip because units do not percolate through the functionality of AIBECS flawlessly at this early stage of development.) We combine the transport operators for DIP and POP into a single tuple via

T_all = (T_DIP, T_POP) ;

(Main.ex-P_model_2_tracers.T_DIP, Main.ex-P_model_2_tracers.T_POP)

Local sources and sinks

Geological Restoring

Because AIBECS will solve for the steady state solution directly without time-stepping the goverining equations to equilibrium, we don't have any opportunity to specify any intial conditions. Initial conditions are how the total amount of conserved elements get specified in most global biogeochemical modelels. Thus to specify the total inventory of P in AIBECS we add a very weak resporing term to, e.g., the DIP equation. The time-scale for this restoring term is chosen to be very long compared to the timescale with which the ocean circulation homogenizes a tracer. Because of this long timescale we call it the geological restoring term, but geochemists who work on geological processes don't like that name! In any event the long timescale allows us to prescribe the total inventory of P in a way that yields the same solution we would have gotten had we time-stepped the model to steady-state with the total inventory prescribed by the initial condition. We define the geological restoring simply as

geores(x, p) = (p.DIPgeo .- x) / p.τgeo

geores (generic function with 1 method)

Uptake (DIP → POP)

For the phosphate uptake, we must take some precautions to avoid the pole of the hyperbolic term. That is, we ensure that uptake only occurs when DIP is positive via the relu function

relu(x) = (x .≥ 0) .* x

relu (generic function with 1 method)

Then we simply write the uptake as a Julia function, following the definition of $U$:

function uptake(DIP, p)
    τ, k, z₀ = p.τ, p.k, p.z₀
    DIP⁺ = relu(DIP)
    return 1/τ * DIP⁺.^2 ./ (DIP⁺ .+ k) .* (z .≤ z₀)
end

uptake (generic function with 1 method)

where we have "unpacked" the parameters to make the code clearer and as close to the mathematical equation as possible.

Remineralization (POP → DIP)

remineralization(POP, p) = p.κ * POP

remineralization (generic function with 1 method)

Net sources and sinks

We lump the sources and sinks into sms functions (for Sources Minus Sinks):

function sms_DIP(DIP, POP, p)
    return -uptake(DIP, p) + remineralization(POP, p) + geores(DIP, p)
end
function sms_POP(DIP, POP, p)
    return  uptake(DIP, p) - remineralization(POP, p)
end

sms_POP (generic function with 1 method)

and create a tuple of those just like for the transport matrices via

sms_all = (sms_DIP, sms_POP) # bundles all the source-sink functions in a tuple

(Main.ex-P_model_2_tracers.sms_DIP, Main.ex-P_model_2_tracers.sms_POP)

For the sake of this notebook, however, we write more efficient BGC functions for the sources and sinks by using "in place" functions that mutate a preallocated vector to avoid allocations. G_DIP! replaces sms_DIP:

function G_DIP!(dDIP, DIP, POP, p)
    τ, k, z₀, κ, DIPgeo, τgeo = p.τ, p.k, p.z₀, p.κ, p.DIPgeo, p.τgeo
    dDIP .= @. -(DIP ≥ 0) / τ * DIP^2 / (DIP + k) * (z ≤ z₀) + κ * POP + (DIPgeo - DIP) / τgeo
    return dDIP
end

G_DIP! (generic function with 1 method)

Note the use of the @. macro to "fuse" the loops together, which is faster than looping for each operation. And G_POP! replaces sms_POP:

function G_POP!(dPOP, DIP, POP, p)
    τ, k, z₀, κ = p.τ, p.k, p.z₀, p.κ
    dPOP .= @. (DIP ≥ 0) / τ * DIP^2 / (DIP + k) * (z ≤ z₀) - κ * POP
    return dPOP
end

G_POP! (generic function with 1 method)

which we bundle into

Gs = (G_DIP!, G_POP!)

(Main.ex-P_model_2_tracers.G_DIP!, Main.ex-P_model_2_tracers.G_POP!)

Parameters

We now define and build the parameters. In AIBECS, this is done by first creating a table of parameters, using the API designed for it. First, we create an empty table:

t = empty_parameter_table()    # empty table of parameters

0 rows × 9 columns

Then, we add every parameter that we need

add_parameter!(t, :DIPgeo, 2.12u"mmol/m^3",
               optimizable = true,
               variance_obs = ustrip(upreferred(0.1 * 2.17u"mmol/m^3"))^2,
               description = "Mean PO₄ concentration")
add_parameter!(t, :τgeo, 1.0u"Myr",
               description = "Geological restoring timescale")
add_parameter!(t, :k, 6.62u"μmol/m^3",
               optimizable = true,
               description = "Half-saturation constant (Michaelis-Menten)")
add_parameter!(t, :z₀, 80.0u"m",
               description = "Depth of the euphotic layer base")
add_parameter!(t, :w₀, 0.64u"m/d",
               optimizable = true,
               description = "Sinking velocity at surface")
add_parameter!(t, :w′, 0.13u"1/d",
               optimizable = true,
               description = "Vertical gradient of sinking velocity")
add_parameter!(t, :κ, 0.19u"1/d",
               optimizable = true,
               description = "Remineralization rate constant (POP to DIP)")
add_parameter!(t, :τ, 236.52u"d",
               optimizable = true,
               description = "Maximum uptake rate timescale")

8 rows × 9 columns

These values of these parameters have been optimized for the same model but embedded in the OCIM1 circulation and should yield a pretty good fit to observations. However, because the circulation has changed from the OCIM1 to the OCIM0.1 used here, the best parameters have probably changed, too. Thus, we will re-optimize these parameters below. But first, we will just try to find the steady-state solution with these initial parameter values.

We generate a new type for the parameters, via

initialize_Parameters_type(t, "Pcycle_Parameters")   # Generate the parameter type

Finally, we can instantiate the parameter vector in a single line:

p = Pcycle_Parameters()

DIPgeo = 2.12e+00 [mmol m⁻³] 
  τgeo = 1.00e+00 [Myr] (fixed)
     k = 6.62e+00 [μmol m⁻³] 
    z₀ = 8.00e+01 [m] (fixed)
    w₀ = 6.40e-01 [m d⁻¹] 
    w′ = 1.30e-01 [d⁻¹] 
     κ = 1.90e-01 [d⁻¹] 
     τ = 2.37e+02 [d] 

where we have used the name of the new type that we created just earlier.

Finding the steady-state solution

Generate state function and Jacobian

We generate the in-place state function F! and its Jacobian in a single line and the out-of-place F for use by the solver. (Using dx = similar(x) creates a vector of the same type and size as x but without spending time copying its values, which is OK because we always fill the entire dx in calls to g_DIP! and G_POP!`.)

nb = length(iwet)
F!, ∇ₓF = inplace_state_function_and_Jacobian(T_all, Gs, nb)
F(x::Vector{Tx}, p::Pcycle_Parameters{Tp}) where {Tx,Tp} = F!(Vector{promote_type(Tx,Tp)}(undef,length(x)),x,p)

F (generic function with 1 method)

Solve for the steady-state

We start from an initial guess,

x = p.DIPgeo * ones(2nb) # initial iterate

382338-element Array{Float64,1}:
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 ⋮                    
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004
 0.0021200000000000004

Create an instance of the steady-state problem

prob = SteadyStateProblem(F, ∇ₓF, x, p)

SteadyStateProblem with uType Array{Float64,1}
u0: [0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212  …  0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212, 0.00212]

and solve it

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

382338-element Array{Float64,1}:
 0.0020475554910165607
 0.002057421172117009 
 0.0017072134842503503
 0.001604394745303311 
 0.0014699843915421508
 0.0014133418232435521
 0.0013149985799484015
 0.0014421910936787152
 0.0013064467964472536
 0.0012054523290650314
 ⋮                    
 4.613921330187689e-9 
 2.672543832284385e-9 
 2.8666492055640978e-9
 4.854081489360668e-9 
 5.1262350724626185e-9
 4.906041130877634e-9 
 5.49232446806064e-9  
 5.310574580913026e-9 
 5.4657917053741616e-9

Plotting

For plotting, we first unpack the state

DIP, POP = state_to_tracers(s, nb, 2)

([0.00204756, 0.00205742, 0.00170721, 0.00160439, 0.00146998, 0.00141334, 0.001315, 0.00144219, 0.00130645, 0.00120545  …  0.00135408, 0.00135561, 0.00132192, 0.00132191, 0.00223223, 0.00223441, 0.00222974, 0.00221617, 0.00222811, 0.00223251], [2.132e-5, 2.14231e-5, 1.77648e-5, 1.66908e-5, 1.52867e-5, 1.4695e-5, 1.36678e-5, 1.49964e-5, 1.35784e-5, 1.25235e-5  …  5.02825e-9, 4.61392e-9, 2.67254e-9, 2.86665e-9, 4.85408e-9, 5.12624e-9, 4.90604e-9, 5.49232e-9, 5.31057e-9, 5.46579e-9])

We will plot the concentration of DIP at a given depth horizon

iz = findfirst(grd.depth .> 2000u"m")
iz, grd.depth[iz]

(17, 2141.005498821681 m)

DIP_3D = rearrange_into_3Darray(DIP, wet3D)
DIP_2D = DIP_3D[:,:,iz] * ustrip(1.0u"mol/m^3" |> u"mmol/m^3")
lat, lon = ustrip.(grd.lat), ustrip.(grd.lon)

([-89.0, -87.0, -85.0, -83.0, -81.0, -79.0, -77.0, -75.0, -73.0, -71.0  …  71.0, 73.0, 75.0, 77.0, 79.0, 81.0, 83.0, 85.0, 87.0, 89.0], [1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0  …  341.0, 343.0, 345.0, 347.0, 349.0, 351.0, 353.0, 355.0, 357.0, 359.0])

Create the Cartopy canvas, make the data cyclic in order for the plot to look good in Cartopy, and plot the filled contour via contourf

ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
clf()
ccrs = pyimport("cartopy.crs")
cfeature = pyimport("cartopy.feature")
ax = subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))
ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#CCCCCC")      # gray land
lon_cyc = [lon; 360+lon[1]]
DIP_2D_cyc = hcat(DIP_2D, DIP_2D[:,1])
plt = contourf(lon_cyc, lat, DIP_2D_cyc, levels=0:0.2:3.6, transform=ccrs.PlateCarree(), zorder=-1)
colorbar(plt, orientation="horizontal");
title("PO₄ at $(string(round(typeof(1u"m"),grd.depth[iz]))) depth using the OCIM0.1 circulation")
gcf()

Optimizing the model parameters

To optimize the model, we will use data from the World Ocean Atlas (WOA). Specifically, we will use the observed PO₄ climatology, which the WOA produces every few years. Thanks to the WorldOceanAtlasTools package combined with some AIBECS functionality, this is done a just a few commands:

using WorldOceanAtlasTools
μDIPobs3D, σ²DIPobs3D = WorldOceanAtlasTools.fit_to_grid(grd, 2018, "phosphate", "annual", "1°", "an")
μDIPobs, σ²DIPobs = μDIPobs3D[iwet], σ²DIPobs3D[iwet]
μx = (μDIPobs, missing)
σ²x = (σ²DIPobs, missing)

([Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf  …  Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf], missing)

We define some hyper parameters for the weight of the mismatch of DIP, POP, and the parameters with observations / their prior estimates.

ωs = [1.0, 0.0] # the weight for the mismatch (weight of POP = 0)
ωp = 1e-4       # the weight for the parameters prior estimates

0.0001

We then generate the objective function and some derivatives using the PO₄ mean and variance

v = ustrip.(vector_of_volumes(wet3D, grd))
f   =   generate_objective(ωs, μx, σ²x, v, ωp, mean_obs(p), variance_obs(p))
∇ₓf = generate_∇ₓobjective(ωs, μx, σ²x, v, ωp, mean_obs(p), variance_obs(p))
∇ₚf = generate_∇ₚobjective(ωs, μx, σ²x, v, ωp, mean_obs(p), variance_obs(p))

(::getfield(AIBECS, Symbol("#∇ₚf#84")){Float64,Array{Float64,1},Array{Float64,1}}) (generic function with 1 method)

For the optimization, we will use the F-1 algorithm from the F1Method package.

using F1Method

For the F-1 algorithm, we must initialize the memory buffer

mem = F1Method.initialize_mem(s, p)

F1Method.Mem([0.00204756, 0.00205742, 0.00170721, 0.00160439, 0.00146998, 0.00141334, 0.001315, 0.00144219, 0.00130645, 0.00120545  …  5.02825e-9, 4.61392e-9, 2.67254e-9, 2.86665e-9, 4.85408e-9, 5.12624e-9, 4.90604e-9, 5.49232e-9, 5.31057e-9, 5.46579e-9], nothing, [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0], nothing)

We then define the objective, gradient, and Hessian using the F1Method package API

objective(p) = F1Method.objective(f, F, ∇ₓF, mem, p, CTKAlg())
gradient(p) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, p, CTKAlg())
hessian(p) = F1Method.hessian(f, F, ∇ₓf, ∇ₓF, mem, p, CTKAlg())

hessian (generic function with 1 method)

Because we want the parameters to remain positive, we define the change of variables The change variable from vector λ = log(p) to Parameters p:

λ2p(λ) = AIBECS.opt_para(exp.(λ))

λ2p (generic function with 1 method)

and the corresponding derivatives:

∇λ2p(λ) = exp.(λ)
∇²λ2p(λ) = exp.(λ)

∇²λ2p (generic function with 1 method)

Reverse change of variables from Parameters p to vector λ

p2λ(p) = log.(optvec(p))
λ = p2λ(p)

obj(λ) = objective(λ2p(λ))
using LinearAlgebra
grad(λ) = gradient(λ2p(λ)) * Diagonal(∇λ2p(λ))
function hess(λ)
    ∇p = Diagonal(∇λ2p(λ)) # for variable change
    ∇²p = Diagonal(∇²λ2p(λ)) # for variable change
    G = vec(gradient(λ2p(λ)))
    H = hessian(λ2p(λ))
    return ∇p * H * ∇p + Diagonal(G) * ∇²p
end

hess (generic function with 1 method)

Use optim to optimize

using Optim

Add storage for use with Optim

m = length(p)
grad(s, λ) = s[1:m] .= vec(grad(λ))
hess(s, λ) = s[1:m,1:m] .= hess(λ)

hess (generic function with 2 methods)

Set options of Optim

opt = Optim.Options(store_trace = false, show_trace = true, extended_trace = false)

                x_abstol = 0.0
                x_reltol = 0.0
                f_abstol = 0.0
                f_reltol = 0.0
                g_abstol = 1.0e-8
                g_reltol = 1.0e-8
          outer_x_abstol = 0.0
          outer_x_reltol = 0.0
          outer_f_abstol = 0.0
          outer_f_reltol = 0.0
          outer_g_abstol = 1.0e-8
          outer_g_reltol = 1.0e-8
           f_calls_limit = 0
           g_calls_limit = 0
           h_calls_limit = 0
       allow_f_increases = false
 allow_outer_f_increases = false
        successive_f_tol = 0
              iterations = 1000
        outer_iterations = 1000
             store_trace = false
           trace_simplex = false
              show_trace = true
          extended_trace = false
              show_every = 1
                callback = nothing
              time_limit = NaN

The starting parameters are

p

DIPgeo = 2.12e+00 [mmol m⁻³] 
  τgeo = 1.00e+00 [Myr] (fixed)
     k = 6.62e+00 [μmol m⁻³] 
    z₀ = 8.00e+01 [m] (fixed)
    w₀ = 6.40e-01 [m d⁻¹] 
    w′ = 1.30e-01 [d⁻¹] 
     κ = 1.90e-01 [d⁻¹] 
     τ = 2.37e+02 [d] 

Run optimization

results = optimize(obj, grad, hess, λ, NewtonTrustRegion(), opt)

Because we started from an other optimized set of parameters, this should only take a few iterations. The optimized parameters are

p_optimized = λ2p(results.minimizer)

After optimizing these parameters, we expect the fit to be better. Let's have a look at the optimized steady-state solution at the same depth as earlier, and compare it to the observed DIP field. First, let's calculate the steady-state with the optimized parameters.

prob_optimized = SteadyStateProblem(F, ∇ₓF, s, p_optimized)
s_optimized = solve(prob_optimized, CTKAlg()).u

Like earlier, we fill in a 3D array of NaNs with the values we want to plot. Here, we fill the arrays with the fractional difference bewteen modeled and observed DIP:

DIPold, _ = state_to_tracers(s, nb, 2)
DIPnew, _ = state_to_tracers(s_optimized, nb, 2)
δDIPold_3D = fill(NaN, size(grd))
δDIPold_3D[iwet] .= 100(DIPold - μDIPobs) ./ μDIPobs
δDIPnew_3D = fill(NaN, size(grd))
δDIPnew_3D[iwet] .= 100(DIPnew - μDIPobs) ./ μDIPobs

And take a slice at depth index iz

δDIPold_2D = δDIPold_3D[:,:,iz]
δDIPnew_2D = δDIPnew_3D[:,:,iz]

We make those 2D slices cyclic along the longitude for Cartopy

δDIPold_cyc = hcat(δDIPold_2D, δDIPold_2D[:,1])
δDIPnew_cyc = hcat(δDIPnew_2D, δDIPnew_2D[:,1])

And we plot the mismatch of the old, OCIM1-optimized field

figure()
δDIPlevels = -20:2:20
ax = subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))
ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#CCCCCC")      # gray land
plt2 = contourf(lon_cyc, lat, δDIPold_cyc, cmap="PiYG_r", levels=δDIPlevels, transform=ccrs.PlateCarree(), zorder=-1, extend="both")
cbar2 = colorbar(plt2, orientation="horizontal", extend="both")
cbar2.set_label("δDIP / DIP [%]")
title("old DIP mismatch at $(string(round(typeof(1u"m"),grd.depth[iz]))) depth using the OCIM0.1 circulation")
gcf()

We can see a significant negative bias in the North Atlantic. This is probably due to having swapped the circulation from OCIM1 to OCIM0.1. Let's plot the OCIM-0.1-optimized mismatch

figure()
ax = subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))
ax.add_feature(cfeature.COASTLINE, edgecolor="#000000") # black coast lines
ax.add_feature(cfeature.LAND, facecolor="#CCCCCC")      # gray land
plt3 = contourf(lon_cyc, lat, δDIPnew_cyc, cmap="PiYG_r", levels=δDIPlevels, transform=ccrs.PlateCarree(), zorder=-1, extend="both")
cbar3 = colorbar(plt3, orientation="horizontal", extend="both")
cbar3.set_label("δDIP / DIP [%]")
title("new DIP mismatch at $(string(round(typeof(1u"m"),grd.depth[iz]))) depth using the OCIM0.1 circulation")
gcf()

This shows that most of the negative bias has been corrected. (The green patch in the North Atlantic has dissipated.) Of course, this optimization was "easy", because there were only 2 tracers, few parameters, and we started from a pretty good solution. However, this type of optimization can really improve the skill of models, even current "state-of-the-art" ones. Such optimizations turn average models into better tools for cutting-edge oceanography research.

This page was generated using Literate.jl.