Structure of Thermodynamics vs Dynamics

[simone-silvestri]

In a sea-ice model, thermodynamics and dynamics are intrinsically connected (tracers are advected by velocities and rheologies need tracers).
Therefore, it makes sense to have them as the same level instead of one inside the other.

I think it would be nice to modify a bit the current structure as to have an overarching SeaIceModel that contains the necessary prognostic variables (also a bit discussed in #19)

• velocities
• tracers
• thickness
• concentration

and

• thermodynamics with all the ruleset for freezing/melting/heat exchange (for example SlabIceThermodynamics) that would be calculated in the compute_tracer_tendencies! function
• dynamics with all the rulesets for momentum stepping/rheologies/solvers (for example ExplicitMomentumSolver and ViscoPlasticRheology) calculated in the step_momentum! function

[simone-silvestri]

Then the tracers to be computed will be required by the thermodynamics argument (a little bit like we do closures now in oceanigans). The velocities would be always 2D fields

[simone-silvestri]

Also, this structure generalizes nicely to extension to, for example, a snow model

[glwagner]

This makes sense. Here's the current structure of SlabSeaIceModel:

struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS}
    grid :: GR
    clock :: CL
    timestepper :: TS
    # State
    ice_thickness :: IT
    ice_concentration :: IC
    top_surface_temperature :: ST
    ice_salinity :: IS
    velocities :: U
    # Boundary conditions
    external_heat_fluxes :: STF
    heat_boundary_conditions :: TBC
    # Internal flux
    internal_heat_flux :: CF
    # Melting and freezing stuff
    phase_transitions :: P
    ice_consolidation_thickness :: MIT
    # Numerics
    advection :: A
end

I'm not so sure about where tracers belong. For example, some tracers may be 2D (not sure which...) but some are certaintly 3D in a multilayer model, like energy, temperature, salinity, porosity.

There are also abstraction issues to sort out regarding the grid for a multilayer model. It could make sense to use a 3D grid for a multilayer model, where "z" would be a non-dimensional coordinate varying from 0, 1. But another possibility is to use a 2D grid and organize the data differently somehow...

[simone-silvestri]

We could always provide a 3D grid which would fix the number of levels in a multilayer model and then force the 2D variables to live in Nothing in the z-direction

[glwagner]

Indeed, there is no challenge with having 2D variables on a 3D grid. We have designed the system to support that.

The hard question is what is clearest for users. In terms of implementation I don't think it matters much what we choose.

I think we want a 3D grid, it's just annoying that the vertical coordinate is called "z". So we will have to clearly document what that means for a sea ice model.