negative snow/water test, beta

[kmdeck]

Purpose

May close CliMA/ClimaCoupler.jl#472

The default parameters of bucket model make it so the beta function recovers the SLIM model (linear beta, switch occurs at 0.75*bucket capacity, and no snow beta).

To-do

TBD: make beta_snow and beta_bucket different functions
TBD: make them functions the user can pass in (for now, the user can control this with parameters of the function)

Content

• Adds several checks of the following form: W = max(W, 0) in functions where being negative would return an incorrect answer and we want to treat the negative number like zero
• adds tests that ensure the tendencies and water fluxes are zero in the case of zero precip, snow melt, but negative snow and water content (evap = sublimation = 0)
• Employs a beta factor for snow. This has the same form as for the bucket, namely:

if x < x_c
    (x/x_c)^p
else
    1
end

This is what we had before for water, with x_c = 0.75*W_f, x = W, and p=1. We've now set p=2 to damp nonlinearly to zero. The branch is used for the bucket because the land evaporates at the potential rate if the bucket is full enough (beta=1), but then drops as it gets drier. The transition happens at x_c [SLIM supplement text Equation 36].

For snow, the sublimation rate should be the potential rate (beta = 1) until the snow vanishes. To achieve this, we use the same function for beta, but use x_c = 0.1 S_c, where S_c is the parameter used to control the albedo of the land surface. For snow levels about S_c, the albedo is the snow albedo, as S drops far below S_c, it becomes the bareground albedo. My only concern here is that the physical situation could be S<<S_c (patchy, thin snow), and the albedo-> bareground albedo, but the vapor flux -> 0, because we treat the snow cover fraction as zero or 1 in all cases but the albedo. In reality, in this situation, the vapor flux should tend towards the bareground vapor flux.

Review checklist

I have:

• followed the codebase contribution guide: https://clima.github.io/ClimateMachine.jl/latest/Contributing/
• followed the style guide: https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/
• followed the documentation policy: https://github.com/CliMA/policies/wiki/Documentation-Policy
• checked that this PR does not duplicate an open PR.

In the Content, I have included

• relevant unit tests, and integration tests,
• appropriate docstrings on all functions, structs, and modules, and included relevant documentation.

---

• I have read and checked the items on the review checklist.

[LenkaNovak]

LGTM, thank you, @kmdeck ! I only had a couple of minor comments.

[LenkaNovak]

Do we want to make 0.75 and 0.1 tunable parameters? It would be good to comment on these choices in the doc string.

[kmdeck]

We can make these tunable params - and p also?

[kmdeck]

Ill add them to the parameter struct

[LenkaNovak]

Should this be eps(FT), so it's consistent with the heaviside function? (also below)

[kmdeck]

I'm not sure. I think both are consistent in a way:

In the heaviside, we are saying that snow < eps(FT) means the snow cover fraction is treated as zero. Here, we are saying that snow being negative means treat it like zero snow. A negative snow would also give zero snow cover fraction, so those are consistent.

If we treat negative snow as eps(FT) as the snow content, the heaviside will return 0 still, since we compare eps(FT) < eps(FT), which is false. so that seems fine too.

Im not sure! if we want zero snow it almost seems like we should return zero snow. using eps(FT) as the heaviside threshold is in case the snow damps to zero but does not exactly reach == 0.

I

[LenkaNovak]

I meant more for the case of 0 < σS < eps(FT), whether it should be treated as 0 or nonzero. It also appears in partition_surface_fluxes(). That said, it looks like the snow_cover_fraction which is only used to scale evaporation, so I think that's ok, as long as the user doesn't get confused how to treat small values of snow. Maybe a comment would do. :)

[LenkaNovak]

But in that case, should the comment on Line 457 (# Equation (3) of the text, snow melt already multipied by snow_cover_fraction) say max(0, σS)? 🤔

[kmdeck]

ah, I think I understand. You're right that we are being inconsistent in the 0 < σS < eps(FT) regime. In principle, the snow_cover_fraction function can be anything (e.g. sigmoid function), and clipping negative values of σS to zero should still be a valid thing to do (they shouldnt depend on each other). I think they do because some of parameterizations depend on σS (e.g. the albedo, the beta function) and some depend on snow cover fraction. partition_surface_fluxes uses both in the same function.

Here is what I did to try and fix this:

• partition_surface_fluxes now works in units of emitting/absorbing area (snow). This is consistent with how we handle fluxes in the documentation and everywhere else - they are fluxes per emitting area and we multiply by emitting area fraction (1-σ) or σ to get fluxes per unit ground area when used in setting tendencies. σ is now applied outside the partition_surface_fluxes now, in Bucket.jl lines 457, 468, 474, and 437.
• The albedo model (linear interpolation between snow and bareground albedo) fundamentally conflicts with a heaviside function for snow cover fraction, so I left this alone for now. The inconsistency remains but it is O(eps) and doesnt affect any conservation. I left a note in the documentation.
• The other inconsistency arises in the beta function, but this function is independent of the snow cover fraction and implemented for numerical reasons only, so I dont know that they need to be consistent. if 0 < σS < eps(FT), this function will be called for snow but it wont be used because it is multiplied by the SCF.

If we want to use the bucket model as a stepping stone for the full land model for testing different fractional coverages within a grid cell, we can make this consistent and use the SCF for the albedo function, at least.

[LenkaNovak]

That seems reasonable to me. Thank you for digging into this, Kat!

[LenkaNovak]

LGTM, @kmdeck , thank you! I only had one more clarification comment.