ACE: Y_fit_max evaluates a different fit than the one in its comment

[lxmota]

Found while reading how porosity enters J2Erosion compared to the new Permafrost cap model. Reporting rather than fixing, because the right answer depends on which form was actually calibrated, and that is Elyce's call.

What the code does

src/LCM/models/J2Erosion_Def.hpp:436-449:

Y_fit_max(T const x, RealType const y)
{
  // UPDATE 3-5-25: changing fit to be the one derived from using purely experimental values,
  // instead of one using *simulated recreations* of the experiments

  // // x = ice saturation;   y = porosity
  // Y_fit = -0.0419 + -0.2972*x + -0.0418*y + 4.7013*x*y

  // -Overall R2: 0.9175
  // -Expt pt R2: -0.0225
  // -Bound pt R2: 0.9688
  // -->Max val = 4.3204 MPa
  return (-0.0419*y  + -0.2972*x + -0.0418*y + 4.7013*x*y) / (4.3204);
}

The comment on line 442 gives the intercept as a bare constant, -0.0419. The code on line 448 multiplies it by the porosity y. So the expression that is evaluated is

(-0.0837*y - 0.2972*x + 4.7013*x*y) / 4.3204

instead of

(-0.0419 - 0.2972*x - 0.0418*y + 4.7013*x*y) / 4.3204

The intercept is dropped and the y coefficient is doubled. This looks like a transcription slip when the 3-5-25 refit was pasted in, not a deliberate reparametrization: the neighbouring E_fit_max (lines 419-429) transcribes its own comment exactly, intercept included.

Size of the discrepancy

Normalized yield multiplier ny, as coded against as commented:

porosity	ice sat	code	comment	difference
0.0	1.0	-0.0688	-0.0785	+0.0097
0.3	0.5	0.1230	0.1162	+0.0068
0.6	1.0	0.5725	0.5686	+0.0039
1.0	1.0	1.0000	1.0000	0.0000

The two forms coincide exactly at y = 1 and diverge by at most 0.0419/4.3204 = 0.0097 of normalized Y elsewhere, so the effect is a small offset rather than a shape change. At the porosity 0.60 used in every ACE deck in the tree it is under 0.5% of the multiplier. That is why nothing has visibly gone wrong, and also why it will not show up on its own.

The decision

Which of the two is the fit that was regressed against the experimental data? If it is the commented form, this is a one-character fix (-0.0419*y becomes -0.0419) plus a rebaseline of the ACE golds. If the coded form is what was actually intended, the comment should be corrected instead so the next reader does not re-open this.

Happy to land either once you say which.

Related, not filed separately

While in the same function: of the three multipliers unit_fit returns, E and Y are floored (lines 577-581, against ACE Residual Elastic Modulus and ACE Soil Yield Strength, then against zero) but K is not. K_fit_min returns a negative multiplier for any porosity below roughly 0.29 at full ice saturation, and for all porosities as saturation goes to zero, so K = hardening_modulus_ * nk reaches the yield surface

ys = SQ23 * (Y + K*eqps_old + sat_mod*(1 - exp(-sat_exp*eqps_old)));

negative and unclamped. That may be intended softening for thawing soil, given the _min in the name and the negative normalizer, but it is asymmetric with how E and Y are handled and worth a deliberate decision. Mentioning it here rather than opening a third issue since it is the same three-line block; say the word if you would rather track it separately.

cc @ebayats

[lxmota]

Related: #123, the other thing that turned up in this same read. Both live in the porosity handling of J2Erosion, both are quiet rather than fatal, and neither affects any deck in the tree today, so they can be triaged together.