Need some changes to avoid negative h2osoi_ice in UpdateState_TopLayerFluxes

[billsacks]

Brief summary of bug

In UpdateState_TopLayerFluxes, there is a somewhat arbitrary epsilon used to determine if h2osoi_ice and h2osoi_liq are close enough to zero that they should be truncated to zero. If the state remains negative after this truncation, we deem this to be a problem. It seems that this tolerance is occasionally exceeded, leading runs to abort. Since this tolerance is somewhat arbitrary, we will loosen it by an order of magnitude.

General bug information

CTSM version you are using: Recent versions

Does this bug cause significantly incorrect results in the model's science? No

Configurations affected: Unknown

Details of bug

For details, see comments in #1253 starting with #1253 (comment)

See also #988

[billsacks]

From discussion with @olyson and @swensosc - We're going to loosen the tolerance further to 1e-10 in the hopes that we'll avoid this problem moving forward. If this problem crops up again, see code in

CTSM/src/biogeophys/SoilFluxesMod.F90

Lines 307 to 359 in 4e126c7

	! Partition the evaporation from snow/soil surface into liquid evaporation,
	! solid evaporation (sublimation), liquid dew, or solid dew. Note that the variables
	! affected here are all related to the snow subgrid patch only because of the use of qflx_ev_snow.
	! In the situations where there are snow layers or there is snow without an explicit snow layer,
	! the partitioned variables will represent the components of snow evaporation
	! (qflx_ev_snow = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
	! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer).
	! In the case of no snow, qflx_ev_snow has already been set equal to qflx_ev_soil (the evaporation
	! from the subgrid soil patch) and the partitioned variables will then represent evaporation from the
	! subgrid soil patch.
	! In the case of urban columns (and lake columns - see LakeHydrologyMod), there are no subgrid
	! patches and qflx_evap_soi is used. qflx_evap_soi = qflx_liqevap_from_top_layer + qflx_solidevap_from_top_layer
	! - qflx_liqdew_to_top_layer - qflx_soliddew_to_top_layer.
	if (.not. lun%urbpoi(l)) then
	if (qflx_ev_snow(p) >= 0._r8) then
	! for evaporation partitioning between liquid evap and ice sublimation,
	! use the ratio of liquid to (liquid+ice) in the top layer to determine split
	if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
	qflx_liqevap_from_top_layer(p) = max(qflx_ev_snow(p)*(h2osoi_liq(c,j)/ &
	(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
	else
	qflx_liqevap_from_top_layer(p) = 0._r8
	end if
	qflx_solidevap_from_top_layer(p) = qflx_ev_snow(p) - qflx_liqevap_from_top_layer(p)
	else
	if (t_grnd(c) < tfrz) then
	qflx_soliddew_to_top_layer(p) = abs(qflx_ev_snow(p))
	else
	qflx_liqdew_to_top_layer(p) = abs(qflx_ev_snow(p))
	end if
	end if
	else ! Urban columns
	if (qflx_evap_soi(p) >= 0._r8) then
	! for evaporation partitioning between liquid evap and ice sublimation,
	! use the ratio of liquid to (liquid+ice) in the top layer to determine split
	if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
	qflx_liqevap_from_top_layer(p) = max(qflx_evap_soi(p)*(h2osoi_liq(c,j)/ &
	(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
	else
	qflx_liqevap_from_top_layer(p) = 0._r8
	end if
	qflx_solidevap_from_top_layer(p) = qflx_evap_soi(p) - qflx_liqevap_from_top_layer(p)
	else
	if (t_grnd(c) < tfrz) then
	qflx_soliddew_to_top_layer(p) = abs(qflx_evap_soi(p))
	else
	qflx_liqdew_to_top_layer(p) = abs(qflx_evap_soi(p))
	end if
	end if
	end if

and see #1282 for some guidance on what might be needed. We do see some potential for issues if there are big differences between h2osoi_liq and h2osoi_ice (in which case the liquid or ice flux might have relatively large errors relative to the liquid or ice state – since relative errors depend on the largest term in the equation), so we may need to do some corrections to ensure that fluxes aren't leading us to exceed either the liq or ice state values, probably when these fluxes are computed in SoilFluxesMod.

[billsacks]

@olyson and @swensosc - I woke up this morning wondering if I could find a relatively easy solution that would get more to the heart of the problem, and spent a bit of time looking into this - though didn't want to spend too long on it (analyzing, developing & testing), so where I currently stand is about where we were yesterday – that, for now, we should just loosen the tolerance. However, I'm inclined to loosen it even more than we suggested - to 1e-9 rather than 1e-10. Are you comfortable with that? I figure that a reduction in state by a factor of 1 billion is effectively reducing the state to 0, so we might as well set it to exactly 0.

That said, I'll document here the thoughts I had this morning in case we want to come back to this - or in case either of you think this is worth pursuing further here.

I started with the assumption that the problem is arising due to big differences in the original h2osoi_liq vs h2osoi_ice. (Before pursuing this further, it would probably be worth testing that assumption by printing the original h2osoi_liq as well as h2osoi_ice in the problem point.) If so, a solution could be to check if the intent is to bring the total state to 0, and if so, adjust each individual flux to accomplish that. Specifically, it looks like a solution could be: when we see that qflx_ev_snow (or qflx_evap_soi for urban) exceed evaporation_limit, set a logical flag on this patch; if this is set, then have a different code path that sets both liqevap and solidevap to take away the entire liquid or solid state (respectively).

I think that this solution would be relatively easy to implement, but then I noticed another possible reason why we might need a looser tolerance for urban: because the application of evaporation_limit that I think matters here is qflx_evap_soi(p) = qflx_evap_soi(p) - frac_sno_eff(c)*(evaporation_demand - evaporation_limit) rather than just qflx_ev_snow(p) = evaporation_limit.

So that realization led me to feel like maybe we just want to go ahead and loosen the tolerance like we discussed yesterday, since I can see multiple possible reasons for the old tolerance to be exceeded and I'm still not sure it's worth the time to dig into this further. But I wanted to run this by you to make sure you're still happy with this plan.

[olyson]

I printed out h2osoi_liq for the problematic point as part of the error message:

563: c, lev_top(c) = 96 0
563: h2osoi_ice_top_orig = 1.999323742483267E-006
563: h2osoi_ice = -2.100641709190665E-018
563: h2osoi_liq_top_orig = 1.807655353197003E-002
563: h2osoi_liq = 0.000000000000000E+000
563: frac_sno_eff = 1.00000000000000
563: qflx_soliddew_to_top_layerdtime = 0.000000000000000E+000
563: qflx_solidevap_from_top_layerdtime = 1.999323742485367E-006

So, h2osoi_liq_top_orig is about four orders of magnitude larger than h2osoi_ice_top_orig.

[swensosc]

I am still okay with the plan to loosen the tolerances. Thanks for your explanation.

…

On Wed, Apr 5, 2023 at 5:23 PM Bill Sacks ***@***.***> wrote: @olyson <https://github.com/olyson> and @swensosc <https://github.com/swensosc> - I woke up this morning wondering if I could find a relatively easy solution that would get more to the heart of the problem, and spent a bit of time looking into this - though didn't want to spend too long on it (analyzing, developing & testing), so where I currently stand is about where we were yesterday – that, for now, we should just loosen the tolerance. *However, I'm inclined to loosen it even more than we suggested - to 1e-9 rather than 1e-10. Are you comfortable with that?* I figure that a reduction in state by a factor of 1 billion is effectively reducing the state to 0, so we might as well set it to exactly 0. That said, I'll document here the thoughts I had this morning in case we want to come back to this - or in case either of you think this is worth pursuing further here. I started with the assumption that the problem is arising due to big differences in the original h2osoi_liq vs h2osoi_ice. (Before pursuing this further, it would probably be worth testing that assumption by printing the original h2osoi_liq as well as h2osoi_ice in the problem point.) If so, a solution could be to check if the intent is to bring the total state to 0, and if so, adjust each individual flux to accomplish that. Specifically, it looks like a solution could be: when we see that qflx_ev_snow (or qflx_evap_soi for urban) exceed evaporation_limit, set a logical flag on this patch; if this is set, then have a different code path that sets both liqevap and solidevap to take away the entire liquid or solid state (respectively). I think that this solution would be relatively easy to implement, but then I noticed another possible reason why we might need a looser tolerance for urban: because the application of evaporation_limit that I think matters here is qflx_evap_soi(p) = qflx_evap_soi(p) - frac_sno_eff(c)*(evaporation_demand - evaporation_limit) rather than just qflx_ev_snow(p) = evaporation_limit. So that realization led me to feel like maybe we just want to go ahead and loosen the tolerance like we discussed yesterday, since I can see multiple possible reasons for the old tolerance to be exceeded and I'm still not sure it's worth the time to dig into this further. But I wanted to run this by you to make sure you're still happy with this plan. — Reply to this email directly, view it on GitHub <#1979 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGRN57BK43LLL7KZF5PMTN3W7X5HZANCNFSM6AAAAAAWQRD3CM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[billsacks]

As I discussed with @olyson yesterday: based on the numbers he provided, I am going to go with the solution I mentioned in my earlier comment:

I started with the assumption that the problem is arising due to big differences in the original h2osoi_liq vs h2osoi_ice. (Before pursuing this further, it would probably be worth testing that assumption by printing the original h2osoi_liq as well as h2osoi_ice in the problem point.) If so, a solution could be to check if the intent is to bring the total state to 0, and if so, adjust each individual flux to accomplish that. Specifically, it looks like a solution could be: when we see that qflx_ev_snow (or qflx_evap_soi for urban) exceed evaporation_limit, set a logical flag on this patch; if this is set, then have a different code path that sets both liqevap and solidevap to take away the entire liquid or solid state (respectively).

because I think this will be a more robust fix that will help prevent this issue from coming up again in the future.