A first stab at a diagnostic surface temperature implementation

[simone-silvestri]

This PR implements two types:

• BulkTemperature
• SkinTemperature

The first type is supposed to represent the default option we have know: the surface temperature is provided and equal to the temperature of "surface" model (either sea ice or ocean). This temperature is fixed and not iterated over.

The second type represents a diagnostic temperature that is diagnosed by a flux balance between the

• top turbulent surfaces (dependent on the surface temperature and iterated by the solver)
• top radiative fluxes (prescribed)
• outgoing longwave (calculated from the surface temperature in the previous iteration)
• internal fluxes: provided from the SkinTemperature, for the moment only a simple diffusive flux with prescribed diffusivity leading to an explicit expression for temperature, but they can become more complicated.

In this second case the temperature is solved for in the flux_balance_temperature(::SkinTemperature, T_ocean, net_external_flux) function

These types are provided to the similarity_theory as a new field named surface_temperature_type. The surface temperature is then stored in a T_surface field which is added to similarity_theory.fields

There is a bit of a pitfall in this implementation: it is difficult to implement temperature-dependent radiative properties, because radiative properties (that depend on i and j and the grid) are passed already computed to the solver.
I would be against passing i, j and the grid to the surface flux solver (it would entangle the solver with oceananigans so it would be difficult to swap it out with a SurfaceFluxes.jl solver), so we could rethink albedo and emissivity computation.

[simone-silvestri]

A very crude first calculation in experiments/single_column_experiments/os_papa_surface_temperature.jl shows the difference between prescribed and diagnostic surface temperature. In the latter case, the internal fluxes are just a very simple diffusive flux with

$$\kappa \frac{T_s - T_o}{\Delta z}$$

where $\kappa = 10^{-2}$, $T_s$ is the incognita, $T_o$ is the temperature of the first ocean cell, and $\Delta z$ is half the upper cell spacing.
The solid lines are the results of DiagnosticSurfaceTemperature and the dashed ones are the respective variables calculated using PrescribedSurfaceTemperature

single_column_profiles.mp4

There is a large difference in the ocean surface temperature (about 1 degree K), but very little in the fluxes.

[simone-silvestri]

However, the surface diffusivity predicted by CATKE in the above video is $O(1)$, much larger than the $\kappa = 10^{-2} m^2s^{-1}$ that I prescribed. If I use $\kappa = 0.5 m^2s^{-1}$ there is even less difference in the surface temperature and barely any in the heat fluxes.

single_column_profiles.mp4

I need to test this on the whole surface to see if there are maybe hot or cold spots that show more difference. As it stands, the net difference in heat flux is less than 1 $Wm^{-2}K^{-1}$. Difficult to assess the impact on a 10-day timescale.

Probably for sea ice where timescales are faster this might matter more.

[simone-silvestri]

There is actually quite a large difference in net heat flux: This is the same computation as in the generate_surface_fluxes.jl but done also with a DiagnosticSurfaceTemperature with a diffusivity of 0.1.
Here I am showing the difference of net heat flux (Q_prescribed - Q_diagnostic) on the left and surface temperature (Ts_prescribed - Ts_diagnostic) on the right.
The surface temperature difference is upward of 2 degrees and the net heat flux has differences upward of 40 Wm2K-1 (which is quite large). In practice, a diagnostic temperature seems to have much more extreme fluxes, i.e heating regions (where solar radiation peaks) are heating much more and cooling regions are cooling much more. There might be a bug in the solver though, because some regions NaN.

[simone-silvestri]

There were some gross bugs in the implementation (including the sign of the fluxes). Correcting the implementation the results look like this (same plot as before)

So in general a diagnostic temperature seems to be mitigating (slightly) the surface fluxes as it is intuitive it should do.
Surprisingly enough, it seems that the solver converges in average more quickly with a Diagnostic surface temperature

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 67 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (4b0d77d) to head (2ea71e3).

Files with missing lines	Patch %	Lines
...sRealmFluxes/similarity_theory_turbulent_fluxes.jl	0.00%	32 Missing ⚠️
...aIceModels/CrossRealmFluxes/surface_temperature.jl	0.00%	20 Missing ⚠️
...Models/CrossRealmFluxes/atmosphere_ocean_fluxes.jl	0.00%	9 Missing ⚠️
...rc/OceanSeaIceModels/CrossRealmFluxes/radiation.jl	0.00%	6 Missing ⚠️

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #278   +/-   ##
=====================================
  Coverage   0.00%   0.00%           
=====================================
  Files         33      34    +1     
  Lines       1876    1910   +34     
=====================================
- Misses      1876    1910   +34     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.