IceOceanModel for coupled ice-ocean simulations

[glwagner]

We've made some progress implementing a simple coupled ice-ocean example over at CliMA/ClimaSeaIce.jl#18.

This work has illustrated that even a simple coupled ice-ocean setup that can cope with ocean freezing and heating (latent heat transfers from ocean to ice), and salinity fluxes during ice freezing and melting, is fairly complicated.

To organize efforts I'm proposing we start now to implement an IceOceanModel abstraction here. This is a significant, complicated effort.

From where I sit here are the main considerations / properties of IceOceanModel:

• (Prescribed) atmosphere representing a potentially time-evolving atmospheric state.
• ice component, which is a Simulation wrapped around a ClimaSeaIce model (right now we only have SlabSeaIceModel).
• ocean component, which is a Simulation wrapped around Oceananigans.HydrostaticFreeSurfaceModel.
• model(s) for ice-atmosphere fluxes
• model(s) for ocean-atmosphere fluxes
• model(s) for ice-ocean fluxes
• time_step!(iom::IceOceanModel, dt)

There are two directions to take with the constructor for IceOceanModel:

1. Require that the component models are built by the user and then passed to IceOceanModel. For this to work, the user will have to correctly build the boundary conditions for both component models in a way that is compatible with a coupled simulation; for example we require arrays for the ocean surface fluxes. The ice-atmosphere fluxes further would require advance knowledge of the prescribed atmospheric state in the current API.
2. Collect dictionaries of keyword arguments in the IceOceanModel constructor, and use these to construct the component models under the hood.

I think the in the long run, mode 1 is preferred because it provides more flexibility. However, it is more challenging from a software design perspective because it requires a complete interface for extracting necessary information about boundary conditions from the component models. There would necessarily be a lot of input validation as well; eg a simple example is that the ocean model was not constructed with AbstractArray / mutable surface boundary conditions on temperature. We would have to throw sensible errors for those cases.

Given the difficulties and complexities of mode 1, I propose that we begin building software around the "mode 2" user interface, and see where that takes us.

There are quite a few additional design considerations that are hard to anticipate at the moment. For example, there is a trade-off between interpretability and flexibility when designing the interfaces for flux specification. Moreover, we probably want convenience features for building consistent flux specifications between the different components: we don't want to use different Stefan-Boltzmann constants to compute radiative emission from ice surfaces and ocean surfaces. We don't want incoming shortwave to heat up ice, but not the ocean. Incoming shortwave has the special complexity that it can penetrate shallow ice layers and heat the ocean below. Momentum fluxes from air to sea depend on ice rheology --- if we use "free-drift" ice rheology (or nothing rheology), then air-sea fluxes are independent of the presence of ice cover.