Ocean state bounds, limiters and how to deal with out-of-bound values

[alicebarthel]

Hi team,

I have had questions on the Eos clamping and have concerns about using the ocean funnel values as a strong constraint on the ocean state or the eos calculation. I want to make sure we discuss this widely.

The original scope was to implement an ocean funnel check. After concerns during the review process (e.g. #191 (comment)), the scope was extended to clamping (i.e. limiting the range of T,S over which we use the polynomial). In its strictest version, this is equivalent to extending the density calculation to a constant density once we are outside the "funnel bounds".
As @cbegeman mentioned, this may mean that a blob of outside-bound state doesn't get appropriately treated (e.g. if it fails to trigger convective mixing).

@katsmith133 and I had also discussed the need for a separate ocn_validate_state() check in Omega to ensure we are in the physical range that our calculations assume. Relevant issue: #297

Therefore, I think it best to expand the scope of this discussion to :

• decide on the hard bounds of the acceptable ocean state.
  These should be bounds that the state is regularly checked against (e.g. ocn_validate_state()) and if exceeded, the model is expected to throw an error and exit.
• determine how to treat values that fall between the strict ocean funnel (narrow state range) and the hard bounds above, if applicable.
  This could include a regularization of the polynomial outside the funnel bounds. E.g. a linear Eos beyond the bounds may be a better approximation than a constant density.
• describe additional checks (online or offline) or error terms we would want implemented as safeguards.

Suggestions:

• The code could include a check on DrhoDt specifically in/around the EOS code to ensure reasonable behavior. At this stage we do not have access to SpecVol/Rho from the previous time step in the Eos code, so we should discuss how best to implement this, most likely in the timestepping(?).
• Once we agreed on the hard bounds above, an offline study could check that the poly75t remains well-behaved between the funnel bounds and the hard bounds. If the polynomial derivatives are not crazy, we could decide to use the polynomial to our hard bounds (i.e. no need for extra regularization).

Extra information could be gathered, e.g.:

• MPAS-O truncation thresholds
• do other model use teos-10 beyond the funnel values? do they use truncation?

I am pinging people whom have expressed interest in these questions, but feel free to join or invite additional people here.
@vanroekel @cbegeman @katsmith133 @philipwjones @sbrus89 @xylar

[alicebarthel]

For context, the original description of the funnel McDougall, 2003.

Here, in the context of the 2003 25-term approximation of the FH95 Gibbs function:
"We are now in a position to determine a formula for ρ as a function of S, θ, and p that is consistent with the FH95 equation of state. We fit the equation of state over the pressure range from 0 to 8000 db (0 to 80 Mpa), but noting that the deep ocean contains neither freshwater nor warm water, we have concentrated the fit on water of environmental relevance. The salinity range used in the fit is 0–40 psu at the surface, but the range is reduced to 30–40 psu at pressures below 5500 db. The minimum salinity used in the fit varies linearly with pressure from 0 to 30 psu between 0 and 5500 db. Similarly, the maximum potential temperature used in the evaluation of the fit is 33°C at 0 db, varying linearly with pressure thereafter down to 12°C at 5500 db. The minimum potential temperature of data that is included in the evaluation of the fit corresponds to the freezing temperature at a pressure of 500 db. That is, for a given salinity, the minimum potential temperature (with a reference pressure of 0 db) was chosen so that if the fluid parcel was moved to a pressure of 500 db, its in situ temperature was the freezing temperature at that salinity and pressure. The rationale for this choice is that icebergs and ice shelves are rarely if ever deeper than 500 db and there are no sources of cooling away from the sea surface (there are only sources of heat, viz. geothermal heating and the dissipation of mechanical energy). This procedure thus defines the coolest potential temperature that can be found anywhere in unfrozen water. In practice we calculated the freezing point of seawater at 500 db from the polynomial appearing in Fofonoff and Millard (1983), the corresponding potential temperatures (with a reference pressure of 0 db) ranging from −0.36°C at 0 psu to −2.6°C at 40 psu.

The above paragraph describes a “funnel” over which we evaluated our fit to the equation of state, and this funnel is shown in Fig. 2. In order to better control the errors in density and particularly the errors in the thermal expansion coefficient, we included a larger range of potential temperature in the fitting procedure but did not include this larger range in the evaluation of the errors of the fitted equation of state. The minimum potential temperature used in the fit was actually 2°C less than the freezing line described above. These potential temperatures would correspond to frozen water at a pressure of approximately 3000 db, which is clearly unphysical since, for example, there is no water of in situ temperature of −4°C observed at this pressure. The maximum potential temperatures that are included in the fit range up to 40°C at 0 db and up to 15°C at 5500 db (and below), varying linearly with pressure between these two pressure limits. We include the extra temperature data in the fit so as to control the well-known problems that can occur when extrapolating polynomials beyond their range of validity, particularly when the equation of state is differentiated with respect to potential temperature to obtain the thermal expansion coefficient."

"Maps of the error in the 25-term fit show that the maximum errors tend to occur on the boundaries of our funnel, as is illustrated in Fig. 4. We suggest that the errors at large potential temperature are less important than those near the freezing line, since it is at low temperatures that deep waters are formed and subsequently spread throughout the world ocean basins. Because of this fact and because the largest errors often occur near the freezing line, a useful measure of error at a given pressure is the typical error at cold temperatures on the S–θ plane of the funnel."

Their Fig 2: The smaller funnel over which the errors were calculated (focusing on the accuracy at low potential temperatures).

[philipwjones]

As mentioned in the original comment, it is important to keep the input values within the valid range because a polynomial fit can exhibit large excursions/oscillation once you step outside the range of validity. And as they note, even getting close to the bounds can increase the error. In MPAS and POP, the earlier J&M polynomial approximations had fixed T, S ranges so it was easier to clamp at the min/max T and S. If the excursions outside the range are small by some measure, it may be less important how we extend. And extending with a linear EOS will need to be careful to match the TEOS value at the boundary so there are no discontinuities.

Regarding checks on ocean state...it is an unfortunate property of our simulations that we can exceed realistic ranges, whether it's forcing biases or numerical over/undershoots. So if you try to constrain the bounds too tightly, you'll abort too often and have to relax those bounds. Checking the range on the low/freezing end is probably more important than the high end (as they also note in their discussion above). Because we probably can't get away with the range as tight as the EOS valid range, we'll still need to clamp the EOS.

[alicebarthel]

Thanks for expanding on that @philipwjones, I understand that simulation states can (temporarily or more systematically) have to live in not-so-realistic states. It's exactly why I want an open conversation about what is a hard vs soft state constraint. The eos clamping is not difficult to implement if we decide to stick to the funnel bounds, but has raised questions. Since we are talking about careful coupling/conservation, I prefer to bring any "limiters" to the surface rather than have an inner subroutine clip values without people being aware of it.

[cbegeman]

These are the max,min input values to the default (JM) EOS in MPAS-Ocean:
https://github.com/E3SM-Project/E3SM/blob/29b14d892aaca9d3c173478c4a320007349d3876/components/mpas-ocean/src/shared/mpas_ocn_equation_of_state_jm.F#L68-L72

   real (kind=RKIND), parameter :: &
      ocnEqStateTmin = -2.0_RKIND, &! valid pot. temp. range
      ocnEqStateTmax = 40.0_RKIND, &
      ocnEqStateSmin =  0.0_RKIND, &! valid salinity, in psu
      ocnEqStateSmax = 42.0_RKIND

[cbegeman]

Let me know if this seems like more of a distraction than a help here. If we need to use clamping as opposed to a linearized EOS between the funnel and allowable bounds at the low temperature end, can we follow the T-S trajectory that is consistent with frazil formation to get the T-S inputs back to the funnel boundary? I imagine this is similar to where we will end up at the end of a time step after frazil formation occurs.

[alicebarthel]

For reference, the gsw_infunnel check from the current GSW-C library:

int
gsw_infunnel(double sa, double ct, double p)
{
    return !(p > 8000 ||
        sa < 0 ||
        sa > 42 ||
        (p < 500 && ct < gsw_ct_freezing(sa, p, 0)) ||
        (p >= 500 && p < 6500 && sa < p * 5e-3 - 2.5) ||
        (p >= 500 && p < 6500 && ct > (31.66666666666667 - p * 3.333333333333334e-3)) ||
        (p >= 500 && ct < gsw_ct_freezing(sa, 500, 0)) ||
        (p >= 6500 && sa < 30) ||
        (p >= 6500 && ct > 10.0)
    );
}

Brainstorming possible edge cases:

• Upper bound on p (8,000 dbar) --> do we have deep trenches or is the topography truncated?

• Lower bound on s --> I think we agree that negative salinity has no physical meaning. Are we happy with 0 being a hard bound or do we need any margin for numerical reason and keep it as a clamping bound?

• Upper bound on s (42) --> this one is clearly exceeded in multiple locations, e.g. the Med and the Red Sea. Using upper bound means we neglect high-density water.

• Upper-Ocean Lower bound on T: Freezing temperature at local (Sa, P) --> do we expect any super-cooled water at p<500? Maybe around frazil or ice runoff etc? @cbegeman

• Mid-Ocean Upper bound in T: decreases linearly with P; maximum from 30C at 500 dbar to 10 C at 6,500 dbar

• Deep Ocean Lower bound on T: Freezing temperature at (Sa, P=500) --> do we expect any super-cooled water at p> 500?

• VDeep Ocean Lower bound on S: Do we expect case of very freshwater (S<=30) at high pressure (P> 6,500 dbar)? I was thinking of subglacial hydrology and grounding line melt but maybe it does not occur at that high a pressure?

• VDeep Ocean Upper bound on T: maximum of 10 C allowed below 6,500 dbar

[cbegeman]

Are we happy with 0 being a hard bound

Yes

do we have deep trenches or is the topography truncated?

I don't believe topography is truncated but I don't expect there to be many places with p > 8000 so I think this is reasonable. I checked a Icos30 mesh with MALI topography and <200 cells in a ~500k cell mesh are below this depth.

do we expect any super-cooled water at p<500?

Yes

do we expect any super-cooled water at p> 500

Yes

Do we expect case of very freshwater (S<=30) at high pressure (P> 6,500 dbar)?

No. The pressures wouldn't get that high.

[philipwjones]

@alicebarthel Regarding your comment on limiters, one of the arguments for just doing EOS clamping is that we're not changing the model state (T,S), so there are no energy or conservation considerations. It's equivalent to just using a different EOS that still gives consistent values for a given T,S. If we actually change the model state directly, then we do need to account for the energy added/removed to make that change.

[philipwjones]

The implementation is correct, but you are comparing apples and oranges (or maybe gala vs delicious). Think of it this way, I can also write rho(T,S,P) = rho_linear(Tlin, Slin, P), but that doesn't mean we've somehow lost/gained mass. Just means we're using a different EOS. And our new EOS happens to agree with TEOS through the funnel range but returns a different value outside that range so our polynomial approximation is different. If you look at the equations, the only thing that matters is that we have some function rho(T,S,P). It doesn't specify what that function should be.

[alicebarthel]

Sorry if I am being dense(!): are you responding to my concern about density differences, ie. what goes into the momentum equation, or my point about mis-accounting?
In the case of the mis-accounting, I was mostly thinking of whether a global sum() would be affected by an existing out-of-bounds region. Now that I think about it, it might be a confusion on my part as to whether we expect to conserve sum(S), sum(hS), or sum(rhohS). In the latter case, I can see how using either EOS may give us consistent results, as long as rho, h, and S are all consistent with each other locally.

Update after in-person discussion: it sounds like the question about density differences still remains.

[philipwjones]

Right, the discussion so far was on the conservation/accounting aspects. But we still need to decide how best to extend the EOS outside the funnel for physical realism and impact on pressure gradient and buoyancy. On the low end, I suspect we're only outside the funnel if somehow we've dropped below the freezing temp for a given P, S (though S maybe also out of range). So if we have frazil ice formation on, we should never exceed that bound as long as we're using the TEOS freezing temp in the frazil calculation. If we don't have frazil on, then I think it still makes sense have the density remain at the freezing point and it's fine that the density doesn't vary much. On the high end, it might make sense to come up with a better way of extending.

[philipwjones]

It might also be worth looking at an existing (MPAS) simulation to see just how often, and by how much, we fall outside the TEOS funnel. Just to get an idea what the pathologies typically are.

[cbegeman]

So if we have frazil ice formation on, we should never exceed that bound as long as we're using the TEOS freezing temp in the frazil calculation.

I just want to clarify that, at least in MPAS-O on a per-time-step basis, we can drop below the freezing temperature. The frazil tendency is based on the freezing point determined using the state at the previous time step. When the tracer tendencies are applied toward the end of the time step, they may bring the temperature below the freezing point. The density is then determined and there is no additional frazil formation (except via coupler corrections) before the pressure gradient term is determined on the next time step (using the state from the previous time step).

[alicebarthel]

Summary of today's discussion:
To balance rigor and modeling pragmatism, and given that other models (MOM6, Oceanigans) use TEOS-10 without hard clamping, we agreed to 1) print warnings when the EOS is used outside of its intended range, 2) have hard exits through an ocean_validate_state().
We can implement stronger limits or clamping if we get a concerning number of warnings.

Action items:

• Put ocean_validate() type testing on the roadmap
• Switch existing clamping to an option
• Add Eos checks that print out warnings for the "cube" and "funnel" ranges

Additional checks could include limiters on individual tendencies, including the drho fron the eos.

[philipwjones]

Sorry I missed the discussion. I'm fine with the approach. Because this is a performance-sensitive bit of code, we'll need to explore the most efficient way to perform the checks without inhibiting performance too much.