Climatemachine.jl: Some temperature profiles lead to imbalanced states

[mwarusz]

Description

Some temperature profile choices lead to imbalanced reference states where the imbalance
is measured as max |dp/dz + rho * g|. For example, I get the following results as a function of the number of vertical elements

Profile = DecayingTemperatureProfile{Float64}(290.0, 290.0, 8484.271021693852)
Ne_z = 2, max |dp/dz + ρ * g| = 1.3063884178876606e-6
Ne_z = 4, max |dp/dz + ρ * g| = 1.666213300666186e-8
Ne_z = 8, max |dp/dz + ρ * g| = 5.654126056242603e-10
Profile = DecayingTemperatureProfile{Float64}(300.0, 300.0, 8484.271021693852)
Ne_z = 2, max |dp/dz + ρ * g| = 0.398088286855824
Ne_z = 4, max |dp/dz + ρ * g| = 0.39808957658210886
Ne_z = 8, max |dp/dz + ρ * g| = 0.39808959267902644
Profile = DecayingTemperatureProfile{Float64}(290.0, 220.0, 8484.271021693852)
Ne_z = 2, max |dp/dz + ρ * g| = 0.8197871304156674
Ne_z = 4, max |dp/dz + ρ * g| = 0.8197933363775505
Ne_z = 8, max |dp/dz + ρ * g| = 0.8197934073661308
Profile = DryAdiabaticProfile{Float64}(290.0, 220.0)
Ne_z = 2, max |dp/dz + ρ * g| = 8.830506104118285e-9
Ne_z = 4, max |dp/dz + ρ * g| = 1.1303313840471674e-10
Ne_z = 8, max |dp/dz + ρ * g| = 5.490718990586174e-12

The first and last profiles show the expected behavior. One clue is that the
scale height doesn't seem to depend on temperature.

Steps to reproduce

Run test/Atmos/check-hydrostatic-balance.jl from the branch mw/check-hydrostatic-balance.

Expected behavior

Imbalance is small and gets smaller with mesh refinement.

Additional context

May be causing https://github.com/CliMA/ClimateMachine.jl/issues/1169

For CLIMA developers

• [ ] The bug can be reproduced
• [ ] Labels have been added and developers assigned
• [ ] The bug can be fixed (if not, please explain why not in a comment below)
• [ ] There is a timeline to fix the bug
• [ ] The bug has been fixed (Please link the PR)

[tapios]

Regarding the scale height, the issue here is that l. 126 of TemperatureProfiles.jl,

H_t::FT = FT(R_d(param_set)) * FT(T_surf_ref(param_set)) / FT(grav(param_set)),

should have _T_virt_surf in place of T_surf_ref. But this should not cause deviations from hydrostatic balance. Hydrostatic balance should be satisfied for any value of the scale height.

[tapios]

Your first and second example should both be isothermal, with slightly different temperatures. Something must be going wrong in the pressure calculation in l. 146:

p = -z - profile.H_t * ΔTv′ * log(cosh(z / profile.H_t) - atanh(ΔTv′)) p /= profile.H_t * (1 - ΔTv′^2)

Could it be that _T_virt_surf in l. 124 reverts to the defaults from the parameters, rather than to the input parameter (and this just works out for 290 K as input because this is equal to the default)?

@charleskawczynski ?

[tapios]

If the above hypothesis is right, you should get wrong results for any surface temperature other than 290 K.

[tapios]

And, @charleskawczynski, we need to add a unit test to the reference profiles.

[charleskawczynski]

Adding these tests in this branch. Perhaps we can continue on these tests this week.

[charleskawczynski]

Your first and second example...

Yes, I think both args are optional, and we should take more care with the args here. Some of this may be fixed in a branch in progress.

[bischtob]

@mwarusz, thanks for tracking this issue. It looks like this can easily cause the issues we have been having. I will test some of this before our meeting today as well!

[bischtob]

I can also confirm that using DecayingTemperatureProfile as initial condition WITHOUT using a reference state in a solid-body rotation test case leads to deviation from balanced flow.

[akshaysridhar]

DryAdiabaticProfile allows conservation of mass, energy to machine precision in unperturbed LESConfiguration problem. DecayingTemperatureProfile does not. Summarised in the gist.

https://gist.github.com/akshaysridhar/96eb8d3366ce71b35b35ecd9b2f368a2

[bischtob]

@mwarusz, @akshaysridhar, @thomasgibson, I haven't tested Charlie's fixes yet (#1186), but I have run the solid-body rotation test with a balance initial condition and balance reference state. The run is here, but does not remain at rest.

[mwarusz]

@bischtob
I ran a few simulations using your script. Since it uses explicit time integration I changed cfl_direction to EveryDirection and set CFL to 1 to have better control over stability. I did a mesh refinement study in the vertical and saw the following:

Ne_vert norm(Q) / norm(Q_0)
3       9.9998916222246226e-01
6       9.9999993616241767e-01
12      9.9999999974946929e-01

I tried applying the filter on perturbations from the initial state (not the reference state) and got

Ne_vert norm(Q) / norm(Q_0)
3       1.0000004583408963e+00
6       1.0000000006969725e+00
12      1.0000000000124525e+00

In general, the problem seems to be that, even with analytically balanced states, the hydrostatic balance is only very roughly satisfied for an atmosphere with 45 km top and low resolution in the vertical. To be more precise, the amount of imbalance is max |dp/dz + g * rho| = 0.004 Pa/m for isothermal atmosphere at 280 K with 3 elements and 45 km top. Compare that to the rising bubble setup with 10 elements and 2.5 km top which has max |dp/dz + g * rho| = 1e-11 Pa/m.

We had some ideas at NPS how to make initial/reference states more balanced in the discrete sense. I am interested in trying them.

[tapios]

A simple fix to get discretely balanced reference states would be to replace the calculation of density in ref_state.jl (from ideal gas law) by a discrete calculation from hydrostatic balance: rho = -1/g (dp/dz). Everything else about the reference state can stay as it is.

[charleskawczynski]

A simple fix to get discretely balanced reference states would be to replace the calculation of density in ref_state.jl (from ideal gas law) by a discrete calculation from hydrostatic balance: rho = -1/g (dp/dz). Everything else about the reference state can stay as it is.

This has worked for me in the past, too.

[bischtob]

Closed and continued in #1247.