Horizontal convection with nonlinear equation of state

[navidcy]

Hey, some people from my group were intrigued with the horizonal convection simulation. I was asked whether it's possible to have a nonlinear equation of state and simulate tracers T and S and then have buoyancy b be given as a function of T, S and pressure (or depth)?

The idea is to start with or impose a temperature-salinity configuration so that the surface buoyancy is homogeneous. Then T and S will start diffuse but as they reach certain depth, nonlinearities in the equation of state might give rise to APE and thus circulation.

I'm cc'ing @kialstewart who probably can explain it better. We don't require the full TEOS-10 for this, just the first few polynomials in the polynomial expansion of TEOS.

[glwagner]

This is definitely possible and might be a good thing to add to #1847 ...

The model in #1847 is built with these lines:

model = IncompressibleModel(
           architecture = CPU(),
                   grid = grid,
              advection = WENO5(),
            timestepper = :RungeKutta3,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
                closure = IsotropicDiffusivity(ν=ν, κ=κ),
    boundary_conditions = (b=b_bcs,)
    )

To use conservative temperature and absolute salinity (with a linear equation of state) we write

tracers = (:T, :S)
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4))

as, for example, the ocean wind mixing and convection example does:

Oceananigans.jl/examples/ocean_wind_mixing_and_convection.jl

Line 72 in 89a6c6c

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4))

Nonlinear equations of state defined by SeawaterPolynomials.jl can also be used with SeawaterBuoyancy. To use TEOS10 for example,

using SeawaterPolynomials: TEOS10EquationOfState
buoyancy = SeawaterBuoyancy(equation_of_state=TEOS10EquationOfState())

We also support a few sets of "second-order" polynomials that were introduced by Roquet et al 2015:

using SeawaterPolynomials: RoquetEquationOfStateationOfState
buoyancy = SeawaterBuoyancy(equation_of_state=RoquetEquationOfState(coefficient_set=:SecondOrder)) 

Docs are limited (please help!) but this docstring should explain:

"""
    RoquetSeawaterPolynomial([FT=Float64,] coefficient_set=:SecondOrder)
Returns a `SecondOrderSeawaterPolynomial` with coefficients optimized by
> Roquet et al., "Defining a Simplified yet 'Realistic' Equation of State for Seawater", 
  Journal of Physical Oceanography (2015).
The `coefficient_set` is a symbol or string that selects one of the "sets" of 
optimized second order coefficients. 
Coefficient sets
================
    - `:Linear`: a linear equation of state, `ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ`
    - `:Cabbeling`: includes quadratic temperature term,
                    `ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2`
    - `:CabbelingThermobaricity`: includes 'thermobaricity' term,
                                   `ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2 - R₀₁₁ * Θ * Z`
    - `:Freezing`: same as `:cabbeling_thermobaricity` with modified constants to increase
                   accuracy near freezing
    - `:SecondOrder`: includes quadratic salinity, halibaricity, and thermohaline term,
                       `ρ = ρᵣ + R₁₀₀ * Θ + R₀₁₀ * Sᴬ + R₀₂₀ * Θ^2 - R₀₁₁ * T * Z`
                             + R₂₀₀ * Sᴬ^2 - R₁₀₁ * Sᴬ * Z + R₁₁₀ * Sᴬ * Θ`
The optimized coefficients are reported in 
Table 3 of Roquet et al., "Defining a Simplified yet 'Realistic' 
Equation of State for Seawater", Journal of Physical Oceanography (2015), and
further discussed around equations (12)--(15). The optimization minimizes
errors in estimated horizontal density gradient estiamted from climatological temperature
and salinity distributions between the 5 simplified forms chosen by Roquet et. al
and the full-fledged [TEOS-10](http://www.teos-10.org) equation of state.
"""

Note that SeawaterBuoyancy also supports setting one of temperature or salinity constant tracer, eg:

using SeawaterPolynomials: RoquetEquationOfStateationOfState

tracers = :T
roquet_eos = RoquetEquationOfState(coefficient_set=:SecondOrder)
buoyancy = SeawaterBuoyancy(equation_of_state=roquet_eos, constant_salinity=35.0) 

There's also a reference density associated with both TEOS10EquationOfState and RoquetEquationOfState which can be set with the kwarg reference_density.

Here's a link to the code:

https://github.com/CliMA/SeawaterPolynomials.jl/blob/75c6c1a3cf8b719841e9a4e12994fd906903e689/src/SecondOrderSeawaterPolynomials.jl#L67

[tomchor]

I had no idea we could use SeawaterPolynomials.jl to define different equations of state. I think it would definitely be helpful to include this in one of the examples is the docs. The ocean mixing one is a good candidate in my opinion.

[glwagner]

I think a validation test would be useful too. We haven't run many problems with nonlinear equations of state.

I think this is a simple and neat example:

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.91.063016

This requires figuring out coefficients for a second-order seawater polynomial, though.

[navidcy]

Great!

@kialstewart would the implementation of EOS that the docstring above describes be sufficient? If so, perhaps try to dig out your setup so we have a vague idea on parameter values of interest and perhaps we can try to sit down and set up a simulation when I'm at ANU next week at some point?

[glwagner]

Also of note is that any "second order" polynomial equation of state can be constructed by setting the coefficients in SeawaterPolynomials.jl's SecondOrderSeawaterPolynomial type.

[kialstewart]

Hi @navidcy

Here are a set of values that worked with MITgcm. These describe a 2048264 (nxnynz, so 2-dimensional in x and z) horizontal convection simulation with imposed surface temperature and salinity distributions. The distribution of the surface temperature is sinusoidal in x with a minimum of 5oC at mid-domain and 15oC at the (periodic) endwalls. The distribution of the surface salinity is approximately sinusoidal with a minimum in the middle and maximum at the endwalls; the exact value of the surface salinity at a given location depends on the imposed surface temperature and reference temperature, such that the imposed salinity and temperature fields compensate to give zero density anomaly imposed at the surface. The imposed surface salinity distribution thus depends on the equation of state.

The implementation of the EOS mentioned above would be great. There would be 5 initial EOS configurations to test:

R_100 (linear, only temperature)
R_100 + R_010 (linear, temperature and salinity)
R_100 + R_010 + R_200 (non-linear, cabbeling only)
R_100 + R_010 + R_101 (non-linear, thermobaricity only)
R_100 + R_010 + R_200 + R_101 (non-linear, cabbeling and thermobaricity)

The idiot-check would be that the second configuration produces zero circulation (only a diffusive transport of heat and salt).

Keen to chat more.

Kial

time step
dt = 20.0

number of grid cells
nx = 2048
ny = 2
nz = 64

horizontal grid spacing
dx = 500.0
dy = 500.0

vertical grid spacing (top and bottom, and change smoothly between)
dz_0 = 10.0
dz_n = 100.0

reference temperature (use as initial)
tref = 10.0

reference salinity (use as initial)
sref = (mean of the imposed surface salinity condition)

boundary conditions
no_slip_sides = false
no_slip_bottom = true
rigid_lid = false
implicit_free_surface = true
periodic_x = true
periodic_y = true

vertical diffusivity: temperature (top and bottom values, change smoothly between)
diffKrNrT_0 = 2e-3
diffKrNrT_n = 1e-5

horizontal diffusivity: temperature
diffKhT = 25.0

vertical diffusivity: salinity (top and bottom values, change smoothly between)
diffKrNrS_0 = 2e-3
diffKrNrS_n = 1e-5

horizontal diffusivity: salinity
diffKhS = 25.0

gravity = 9.81

relaxation timescale: temperature
tau_theta_clim_relax = 20

relaxation timescale: salinity
tau_salt_clim_relax = 20