Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Small number protections to CCH water transfer functions #1164

Merged
merged 15 commits into from
May 1, 2024
Merged

Small number protections to CCH water transfer functions #1164

Changes from 2 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
2e8608c
Added protections to Campbell Clapp-Hornberger water transfer functio…
rgknox Feb 22, 2024
f8dc237
reverting the cap on the smoothed cch function, it already has a prot…
rgknox Feb 22, 2024
271282d
Merge branch 'main' into cch-safemath
rgknox Mar 13, 2024
fce4287
Updated error reporting in 1d taylor solve for hydro
rgknox Mar 13, 2024
13c2caa
correct minor build errors
glemieux Mar 17, 2024
8f85de8
Adding ftc protections (ie 0 flux) for perilously low wcs
rgknox Mar 25, 2024
1abba2f
hydro safemath, looking into pinning zero ftc to min psi
rgknox Mar 26, 2024
35b396d
added hydro code to force zero conductance a psi_min
rgknox Mar 29, 2024
57b10d3
Getting zero flux at psi_min in fates hydro
rgknox Apr 2, 2024
a269473
adjusted ordering of initializing bounds on cch functions to avoid un…
rgknox Apr 3, 2024
3a2b877
Merge tag 'sci.1.72.6_api.34.0.0' into cch-safemath-zeroflux
rgknox Apr 5, 2024
cf05475
Added btran to hydro restart, and protecting downscaling rhiz layer f…
rgknox Apr 11, 2024
25fb73a
Updated some comments in low water potential protections
rgknox Apr 29, 2024
a5d9fd9
update to comments in hydro flux disaggregation code
rgknox Apr 29, 2024
ec27207
disaggregating fluxes based on root length only when sumweights is ze…
rgknox May 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 27 additions & 13 deletions biogeophys/FatesHydroWTFMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ module FatesHydroWTFMod
! elastic-caviation region


real(r8), parameter :: min_psi_cch = -9._r8 ! Minimum psi we are willing to track in cch
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you put the unit here (Mpa)? Why 9 instead of 15?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this was actually updated, seems like a github glitch


! Generic class that can be extended to describe
! specific water retention functions

Expand Down Expand Up @@ -710,10 +712,11 @@ subroutine set_wrf_param_cch(this,params_in)
this%th_max = max_sf_interp*this%th_sat
this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max))
this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max))
this%th_min = fates_unset_r8
this%psi_min = fates_unset_r8
this%dpsidth_min = fates_unset_r8

this%psi_min = min_psi_cch
this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max))
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what's the purpose of the tiny(this$th_max)? Wondering if this is an error, since the input to this function should be only water potential.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

originally, this tiny() function was used to prevent a problem with order of operations. Inside these functions we have "if statements" that compare theta and psi to these very psi_min and th_min thresholds we are calculating. By adding or subtracting a tiny, we were preventing the intering of logical blocks that were not yet initialized. But with these changes, we actually don't really use these tiny's as much

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Understand, but does seems to be a little confusing. Maybe we can clean up the codes in later pull requests.

this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_max))

return
end subroutine set_wrf_param_cch

Expand Down Expand Up @@ -751,8 +754,12 @@ function th_from_psi_cch(this,psi) result(th)
if(psi>this%psi_max) then
! Linear range for extreme values
th = this%th_max + (psi-this%psi_max)/this%dpsidth_max
else
th = this%th_sat*(psi/this%psi_sat)**(-1.0_r8/this%beta)
else
if(psi<this%psi_min) then
th = this%th_min
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't this artificially add water to soils?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, but the intention is to go from a value of zero to an inconsequentially small non-zero value. @ckoven and I were talking, this could potentially be problematic if we artificially add small amounts of water and that enables or drives non-zero fluxes that would had otherwise been zero. One possibility would be to force fractional conductivity to zero where it hits theta min. That seems messy though. Honestly, I think if we just use a suitably small minimum value, one that is an incredibly small water content, yet can still accomodate the math needed to switch between psi and theta, this should be sufficient. Perhaps we can try minimum capping at -15 MPa or smaller, perhaps -9 is too high.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will this trigger mass balance check error if we run the model for a long time (let's say 1000 years for spinup) ? Should we should keep track of this error and used it as part of the mass balance error?

Copy link
Contributor Author

@rgknox rgknox Apr 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't believe this is a mass conservation liability. We use these functions to drive fluxes, they never actually over-write the original masses of each pool. It is the fluxes that update the pools. What will happen is that water contents below the psi_min or theta_min values will just end up having zero flux out of that pool, and thus remain unchanged.

For instance here is where we update the actual soil water content, it is based off of fluxes:
https://github.com/NGEET/fates/blob/sci.1.74.0_api.35.0.0/biogeophys/FatesPlantHydraulicsMod.F90#L2729-L2730

else
th = this%th_sat*(psi/this%psi_sat)**(-1.0_r8/this%beta)
end if
end if
return
end function th_from_psi_cch
Expand All @@ -767,8 +774,12 @@ function psi_from_th_cch(this,th) result(psi)

if(th>this%th_max) then
psi = this%psi_max + this%dpsidth_max*(th-max_sf_interp*this%th_sat)
else
psi = this%psi_sat*(th/this%th_sat)**(-this%beta)
else
if(th<this%th_min) then
psi = this%psi_min
else
psi = this%psi_sat*(th/this%th_sat)**(-this%beta)
end if
end if

end function psi_from_th_cch
Expand All @@ -783,9 +794,13 @@ function dpsidth_from_th_cch(this,th) result(dpsidth)

! Differentiate:
if(th>this%th_max) then
dpsidth = this%dpsidth_max
dpsidth = this%dpsidth_max
else
dpsidth = -this%beta*this%psi_sat/this%th_sat * (th/this%th_sat)**(-this%beta-1._r8)
if(th<this%th_min) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good to check on the edge cases

dpsidth = this%dpsidth_min
else
dpsidth = -this%beta*this%psi_sat/this%th_sat * (th/this%th_sat)**(-this%beta-1._r8)
end if
end if

end function dpsidth_from_th_cch
Expand Down Expand Up @@ -915,7 +930,7 @@ subroutine set_wrf_param_smooth_cch(this,params_in)
this%th_max = max_sf_interp*this%th_sat
this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max))
this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max))
this%th_min = 1.e-8_r8
this%th_min = fates_unset_r8
this%psi_min = fates_unset_r8
this%dpsidth_min = fates_unset_r8

Expand Down Expand Up @@ -1037,8 +1052,7 @@ function th_from_psi_smooth_cch(this,psi) result(th)
sat = 1.d0
endif
th = sat * this%th_sat



return
end function th_from_psi_smooth_cch

Expand Down Expand Up @@ -1195,7 +1209,7 @@ function dpsidth_from_th_smooth_cch(this,th) result(dpsidth)
sat_res = 0._r8
alpha = -1._r8/this%psi_sat
lambda = 1._r8/this%beta

pc = 1._r8 * this%psi_from_th(th)
if( pc <= this%scch_pu ) then
! Unsaturated full Brooks-Corey regime.
Expand Down