hydro crashes in sites with no cohorts -> hardening discussion

[rgknox]

When there are no cohorts in a site/column, the hydro code generates nans and other spurious boundary conditions and crashes simulations. As far as we can tell, this first happens in the attempts to create new rhizosphere shell volumes, where the math doesn't work when there are no fine roots at all.

For instance here we calculate total fine root length by looping over cohorts:
https://github.com/NGEET/fates/blob/master/biogeophys/FatesPlantHydraulicsMod.F90#L1714-L1724

If the variable csite_hydr%l_aroot_layer(:) is zero in any layer, there are problems here (local l_aroot is zero and raised to neg power):
https://github.com/NGEET/fates/blob/master/biogeophys/FatesPlantHydraulicsMod.F90#L4207

The solution should be fairly straight forward, as hydro does not need to be run when there are no roots and can have some nominal place-holder values as well as 0 root uptake fluxes.

@pnlfang @JunyanDing have both confirmed this behavior

[mariuslam]

Hi @rgknox ,

Is this issue still actual ?
If yes, I would like to solve this as I think it is relevant for me. I have more than 60% of biomass loss in a year and a 2D solver issue which I think are related. When you say no cohorts at a column level do you mean patch level?

If you think it might be related to my problem I would appreciate more details as I am not sure to understand where changes should be made to overwrite fluxes and not run hydro if no roots.
(@rosiealice, @JunyanDing , @pnlfang )

Thanks a lot in advance :) Marius

[rgknox]

The specific issue brought up in this thread has been addressed. When calculating the shell volumes, we have provisions in place to generate some nominal values when there are no roots present, see here:
https://github.com/NGEET/fates/blob/master/biogeophys/FatesPlantHydraulicsMod.F90#L4412
Do you have any test results that show a crash when there are no cohorts on a column? (as per your question about patch/column, yes its all about how many cohorts on a column because that is the scale within which we resolve soil moisture in CLM/ELM).

[mariuslam]

OK Thanks, Well I don't think that is the problem in the end. I just think the hydro solver struggles when the soil water potential goes lower than -1000MPa.
Is there an easy fix ?

[mariuslam]

With hardening I keep water potential in plant tissues at higher values then default plant compartments and soil potentials

[JunyanDing]

@mariuslam, how do you enforce hardening?
I think too negative the water potential is the cause of your problem. Did the issue go away when you apply hardening?

[mariuslam]

@JunyanDing, (To better show the problem I updated the plot for the soil potential). So the low values for water potential are mostly due to freezing of the 4cm soil layer which reaches -40°C (Both in the reference and in the hardening simulations). In my hardening simulations I decrease Kmax between plant compartments incuding between roots and soil. So I keep the water potential of the plant compartments at values above -15Mpa. I increase the gradient between soil and plant which simply makes the problem occur earlier.
I am running over 30 year atmospheric data, so there is also a link to the biomass building up and influencing soil water potential in the end. because the model crashes the second time it goes through the 6th atmospheric data

[rgknox]

@mariuslam are the units in that plot in MPa? If so, we may have to look through our various functions and perform math on those pressures. It could be that the math operations simply can't handle numbers so large. For instance, in Van Genuchten, we would be raising 150 to a power, that could be bad if the exponents are large. Maybe its not an issue, but worth investigating...

I'm looking through FatesHydroWTFMod.F90 and see some functions where psi is raised to a power that could cause problems (at least with VG). What water retention and conductivity hypothesis did you use?

https://github.com/NGEET/fates/blob/master/biogeophys/FatesHydroWTFMod.F90

Maybe we should identify a safe range of pressures that these functions can operate on, and then get some code written that will protect them somehow. One option would be to have a hard cap on the lower bound of psi. Or.. maybe the problem is something else.

[mariuslam]

Yes Mpa from unchanged output. I reproduced a default run at yakutsk (siberia) with a very recent version : sci.1.48.0_api.17, Showing water potentials reaching -150 in plants and -700 in soil. Potentially much lower when it crashes. I use TFS for plants and CCH for soil. The hardening scheme can keep water potentials high in plants, but would be amazing to have a limit on the soil psi

[ckoven]

I wonder if there is a way to stabilize this by putting a floor on the amount of supercooled liquid water allowed to exist in the soils when they get far below the freezign point? Maybe by putting a max(t_soisno(c,j), min_temp_soilfreeze) instead of the two instances of t_soisno(c,j) here: https://github.com/ESCOMP/CTSM/blob/master/src/biogeophys/SoilTemperatureMod.F90#L1218

with a value of 263K or something like that for min_temp_soilfreeze to let the soil matric potential get low but not super-super low?

[mariuslam]

Thanks Charlie, the solver liked the suggestion :) Here some plots to illustrate what the correction did if I put min_temp_soilfreeze = 263. I am wondering if this value should be a little lower (258K), and if there shouldn't be some complementary fix? It seems like psi still goes to -4000MPa during that extreme year.

[mariuslam]

@rgknox , is there something here (#737) that would solve the low water potentials in the soil ?

[rgknox]

I don't think so, those changes are all really addressing Van Genuchten, which you don't use.

Functionally, -8000 MPa is an amazingly small amount of water. I'm imagining someone trying to reproduce that in a lab, it seems like one atom (edit: molecule, sorry) is alone in that cubic meter of soil, refusing to leave, making its last stand, breaking our math. Even -20MPa is bone dry, nothing there...

I'm thinking a few things. First, we either want to look at the conductance and P-V functions, and make sure they can operate at very low PSI values. At first, I was thinking that we should enforce a cap on psi, and not let it go below. But I think this would cause problems, because you will not get the same volume back if you translate V -> P, and then P->V. Ie capped functions are not bi-directional. I suppose we could put in a linear regime below -20MPa, and give it some artificial weak slope that approximates a cap but allows the functions to be reversable.

Or, it could be an issue with the FTC. On the plant side, we do have a minimum FTC. It was set really low to avoid divide by zeros. But maybe we need to set it even lower. I could be that this minimum fraction of total conductance is slowly generating the steep gradients and forcing water fluxes from pools that are incredibly dry... Maybe try setting it to something incredibly small, like 1e-20?

See here:
https://github.com/NGEET/fates/blob/master/biogeophys/FatesHydroWTFMod.F90#L31

The FTC functions themselves could also be susceptable to underflow at very very large negative PSI... I'm not sure what a good psi_sat value is, but the equation for the derivative has the following term in there (approximating exponent):

dftc_dpsi = a * (-8000 / psi_sat)^-3.5

And playing around with values of psi_sat (which should be negative small), I'm easily getting ftc values << 1e-25.

[mariuslam]

Hi @rgknox,

Thanks for the tips! I tried the easiest first: I set the minimum ftc at 1e-20 and saw no difference, so I played a bit with psi_sat.
Please have a look at the plots below. It does reduce the amplitude of the low water potential spikes, but it does it with all of them. I am not sure whether that is what we want right? The first soil layer temperature sure varies between years, some years it reaches -40, some -20 etc. We should keep I think a strong inter-annual variability, but not as big...
These runs also include the fix proposed by @ckoven. Maybe I will get what I want by downplaying or increase that fix in favor of this one, or do we have no choice but to have a linear interpolation?

[mariuslam]

I have played with the variables mentioned by @rgknox and @ckoven above.
If we want to reduce the extreme peaks I thing there might be no other solution than a linear interpolation in the PV curve. I have no Idea yet of how to proceed in doing that so all tips are welcome.

[mariuslam]

I made a test with linear interpolation below -25MPa, the model runs and reacts as expected. Does the method I use seems reasonable to you ? @rgknox , @rosiealice .
https://github.com/mariuslam/Fates_single_site/blob/cd6deb9ee5e6258852a01303de5b014d048775aa/biogeophys/FatesHydroWTFMod.F90#L700-L758

function psi_from_th_cch(this,th) result(psi)

class(wrf_type_cch)  :: this
real(r8),intent(in)  :: th
real(r8)             :: psi
real(r8)             :: thmin
thmin=this%th_sat*(-25._r8/this%psi_sat)**(-1.0_r8/this%beta)
if(th>this%th_max) then
    psi = this%psi_max + this%dpsidth_max*(th-max_sf_interp*this%th_sat)
else if (th < thmin)   then      ! marius
    psi = -25._r8 - 10._r8 * (thmin-th)/thmin
    write(fates_log(),*)'check1','psi',psi,'th',th,'thmin',thmin
else
    psi = this%psi_sat*(th/this%th_sat)**(-this%beta)
end if

end function psi_from_th_cch