Should TEOS-10 density anomaly omit the (z-dependent) reference profile?

[glwagner]

Curious whether this is a more appropriate definition of the reference density. In principle, this should have no dynamical effect. But it could affect numerics (maybe).

@xkykai @navidcy @jbisits @simone-silvestri

[codecov[bot]]

Codecov Report

Attention: 1 lines in your changes are missing coverage. Please review.

Comparison is base (6f70669) 81.15% compared to head (9543d7a) 81.15%. Report is 14 commits behind head on main.

Files	Patch %	Lines
src/TEOS10.jl	0.00%	1 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #32 +/- ## ======================================= Coverage 81.15% 81.15% ======================================= Files 3 3 Lines 69 69 ======================================= Hits 56 56 Misses 13 13 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[xkykai]

I think this is the right approach, we should omit the inert effects of compressibility since it has no dynamical importance.

[glwagner]

@jbisits if we make this change, then we also have to change the seawater_density utility that you implemented in Oceananigans. Can you weigh in?

[jbisits]

Sorry for being slow on this!

This change looks good. When calculating the in-situ density in seawater_density I think we would just need to add the reference profile back in. We could add a the reference profile to reference_density for a TEOS10EquationOfState

reference_density(eos::TEOS10EquationOfState) = eos.reference_density + r₀(ζ)

when creating the container for a TEOS10EquationOfState? Though it could be a bit misleading when the user entered reference_density does not match the reference_density in the container.

Or we could add the reference profile to a TEOS10 struct

struct TEOS10EquationOfState{P, FT} <: BoussinesqEquationOfState
    seawater_polynomial :: P
      reference_density :: FT
                     r₀ :: FT

    function TEOS10EquationOfState(seawater_polynomial, reference_density)
        FT = eltype(seawater_polynomial)
        P = typeof(seawater_polynomial)
        reference_density = convert(FT, reference_density)
        r₀ = r₀(ζ)
        return new{P, FT}(seawater_polynomial, reference_density, r₀)
    end
end

Then change

https://github.com/CliMA/SeawaterPolynomials.jl/blob/c8198d3003b9e4987b83789c1dc14a83e3ce41a2/src/SeawaterPolynomials.jl#L80C1-L80C69

@inline ρ(Θ, Sᴬ, Z, eos::TEOS10EquationOfState) = eos.reference_density + eos.r₀ + ρ′(Θ, Sᴬ, Z, eos)

Then seawater_density in Oceananigans.jl should still return the correct density variable for a given BoussinesqEquationOfState.

[glwagner]

Sorry for being slow on this!

This change looks good. When calculating the in-situ density in seawater_density I think we would just need to add the reference profile back in. We could add a the reference profile to reference_density for a TEOS10EquationOfState

reference_density(eos::TEOS10EquationOfState) = eos.reference_density + r₀(ζ)

when creating the container for a TEOS10EquationOfState? Though it could be a bit misleading when the user entered reference_density does not match the reference_density in the container.

Or we could add the reference profile to a TEOS10 struct

struct TEOS10EquationOfState{P, FT} <: BoussinesqEquationOfState
    seawater_polynomial :: P
      reference_density :: FT
                     r₀ :: FT

    function TEOS10EquationOfState(seawater_polynomial, reference_density)
        FT = eltype(seawater_polynomial)
        P = typeof(seawater_polynomial)
        reference_density = convert(FT, reference_density)
        r₀ = r₀(ζ)
        return new{P, FT}(seawater_polynomial, reference_density, r₀)
    end
end

Then change

https://github.com/CliMA/SeawaterPolynomials.jl/blob/c8198d3003b9e4987b83789c1dc14a83e3ce41a2/src/SeawaterPolynomials.jl#L80C1-L80C69

@inline ρ(Θ, Sᴬ, Z, eos::TEOS10EquationOfState) = eos.reference_density + eos.r₀ + ρ′(Θ, Sᴬ, Z, eos)

Then seawater_density in Oceananigans.jl should still return the correct density variable for a given BoussinesqEquationOfState.

The reference profile is not a constant though, it is z-dependent.

I just realized that seawater_density uses the function ρ, which is not changed by this PR:

https://github.com/CliMA/Oceananigans.jl/blob/76c32f65bbf1035686515fcc039a45cafcd8c3ff/src/Models/seawater_density.jl#L10

So, no change is needed in fact?

[jbisits]

The reference profile is not a constant though, it is z-dependent.

Yep so if is was stored in a container it would be a Vector.

I just realized that seawater_density uses the function ρ, which is not changed by this PR

Yes that does look to be the case!

All looks good then.

[glwagner]

@jbisits I think it makes sense to keep the TEOS10 implementation of the compressible part of the reference density profile, since this object is called TEOS10EquationOfState after all.

[glwagner]

@simone-silvestri need approval