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.

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

Conversation

rgknox
Copy link
Contributor

@rgknox rgknox commented Feb 22, 2024

A cap has been created on the smallest suction that we are willing to simulate with CCH (Campbell Clapp-Hornberger) water transfer functions. These methods determine the equivalent water content, so that if values below this water content are passed to the functions, it operates at the lowest threshold instead of breaking. The current threshold is -15 MPa. This is incredibly small, and should generate conductances that indistinguishable from conductances generated with near zero water contents.

A further modification was made in a second round of development here. All water retention curves are now designed to have a minimum value. For CCH it is -15 MPa, for Van Genucten it is the residual water content (existing parameter). We now force the fraction of maximum conductivity to be zero at the minimum value, this currently includes stomatal conductance. This is applied by determining a weighting function that blends in a zero value. The weighting is 1 at the minimum suction, and decreases to zero on an exponential. I used an exponent that attenuates to a value of about 0.175 after increasing 1 MPa above minimum, and 0.025 after 2 MPa above minimum. This function makes the derivative of conductance with suction continuous and smooth, but most importantly, it prevents the flow of water in plant and soil media that are completely parched.

This method of ensuring fluxes are zero at these super low water contents will not affect re-hydration, because conductance between nodes is driven by the suction of the up-stream compartment (ie the one with more total potential, or less suction). However, this will prevent completely desicated compartments from going into negative water content.

Description:

Collaborators:

@olyson @jennykowalcz

Expectation of Answer Changes:

This should prevent crashes, particularly in soil water calls using frozen soils.

Fixes: #1168
Fixes: #1163

Checklist

If this is your first time contributing, please read the CONTRIBUTING document.

All checklist items must be checked to enable merging this pull request:

Contributor

  • The in-code documentation has been updated with descriptive comments
  • The documentation has been assessed to determine if updates are necessary

Integrator

  • FATES PASS/FAIL regression tests were run
  • Evaluation of test results for answer changes was performed and results provided

Documentation

Test Results:

CTSM (or) E3SM (specify which) test hash-tag:

CTSM (or) E3SM (specify which) baseline hash-tag:

FATES baseline hash-tag:

Test Output:

@glemieux glemieux self-requested a review February 26, 2024 18:01
@rgknox
Copy link
Contributor Author

rgknox commented Feb 26, 2024

This may address: #1154

@jennykowalcz
Copy link

This may address: #1154

I tried this out -- in the control case, without this PR, the case fails with this error:

128: PIO: FATAL ERROR: Aborting... An error occured, Writing variables (number of variables = 425) to file (./hydro_nocomp
_2pft_campbell_control.elm.h0.0002-11.nc, ncid=141) using PIO_IOTYPE_PNETCDF iotype failed. Non blocking write for variabl
e (FATES_ROOTWGT_SOILMATPOT, varid=274) failed (Number of subarray requests/regions=1, Size of data local to this process 
= 288). NetCDF: Numeric conversion not representable (err=-60). Aborting since the error handler was set to PIO_INTERNAL_E
RROR... (/global/u2/j/jkowalcz/E3SM-test/E3SM/externals/scorpio/src/clib/pio_darray_int.c: 395)

In the test case, pulling in the changes to FatesHydroWTFMod.F90, it fails with this error:

 27:  Site plant water balance does not close
 27:  delta plant storage:  -4.101476403268849E-005  [kg/m2]
 27:  integrated root flux:  -2.537539597419707E-002  [kg/m2]
 27:  transpiration flux:   3.019638043594088E-005  [kg/m2]
 27:  end storage:   1.310177109332794E-002
 27:  pre_h2oveg  1.314278585736063E-002
 27:  ENDRUN:
 27:  ERROR in FatesPlantHydraulicsMod.F90 at line 2808  

In both cases I'm using this E3SM branch with this FATES branch (with the change to FatesHydroWTFMod for the test case). I can try also testing with this PR and the hydro stability PR.

The region used in these runs is a lat/lon rectangle from 15 N - 25 S and 30 W - 85 W so it includes very dry regions outside the Amazon basin where FATES-Hydro is unhappy.

@glemieux
Copy link
Contributor

glemieux commented Mar 6, 2024

Running regression tests on derecho I found that while this resolves #1163, I'm seeing a run failure due to the 1D hydro solver failing with a large wb_step_error:

2461  Could not find a stable solution for hydro 1D solve
2462 
2463  error code:            1
2464  error diag:    0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000
2465  lat:   -7.0000000000000000      longitidue:   305.00000000000000
2466  is recruitment:  F
2467  layer:            7
2468  wb_step_err =    1.6512024553098321E+139
2469  q_top_eff*dt_step =    0.0000000000000000
2470  w_tot_beg =    102.64162586485499
2471  w_tot_end =    1.6512024553098321E+139
2472  leaf water:    6.3082104539614020E-002  kg/plant
2473  stem_water:    5.2044616695225525E-002  kg/plant
2474  troot_water:    1.6178663624130458E-002
2475  aroot_water:    1.8182979624968509E-002
2476  LWP:   -2.1266161551233220
2477  dbh:   0.40361081540638905
2478  pft:            6
2479  z nodes:    1.2616446404957407       0.60582232024787031      -0.19099388737231493      -0.79999999999999993      -0.79999999999999993
2480  psi_z:    1.2364117476858283E-002   5.9370587384291323E-003  -1.8717400962486863E-003  -7.8399999999999997E-003   0.0000000000000000
2481  vol,    theta,   H,  Psi,     kmax-
2482  flux:             0.0000000000000000
2483  l:   8.1613119481785556E-005  0.77294073477601488        3.3911760124701393        3.3788118949932811
2484                           4.3843160661903065E-004
2485  s:   8.0363738039128391E-005  0.64761318929546885       -3.6111596928010109E-002  -4.2048655666439241E-002
2486                           1.2058111921054136E-003
2487  t:   2.4966046503811633E-005  0.64802665578870955       -3.6564678574268254E-002
2488                           8.0381337022338061E-004
2489  a:   1.3824320409604611E-006  0.75035604020636570       -2.2154089112156015E-003
2490                     in:   1.6033112452206815E-003
2491                    out:   3.0709050878900186E-006
2492  r1:  0.21880148501536775       0.46850359885155290        6.9272904988199870E+156

I'm guessing that this is due to something causing dth_node to blow up within Hydraulics_Tridiagonal?

Note that the log stops writing early due to #1168, but I'm not entirely sure if that error is actually due to the above error.

@jennykowalcz
Copy link

This may address: #1154

I tried this out -- in the control case, without this PR, the case fails with this error:

128: PIO: FATAL ERROR: Aborting... An error occured, Writing variables (number of variables = 425) to file (./hydro_nocomp
_2pft_campbell_control.elm.h0.0002-11.nc, ncid=141) using PIO_IOTYPE_PNETCDF iotype failed. Non blocking write for variabl
e (FATES_ROOTWGT_SOILMATPOT, varid=274) failed (Number of subarray requests/regions=1, Size of data local to this process 
= 288). NetCDF: Numeric conversion not representable (err=-60). Aborting since the error handler was set to PIO_INTERNAL_E
RROR... (/global/u2/j/jkowalcz/E3SM-test/E3SM/externals/scorpio/src/clib/pio_darray_int.c: 395)

In the test case, pulling in the changes to FatesHydroWTFMod.F90, it fails with this error:

 27:  Site plant water balance does not close
 27:  delta plant storage:  -4.101476403268849E-005  [kg/m2]
 27:  integrated root flux:  -2.537539597419707E-002  [kg/m2]
 27:  transpiration flux:   3.019638043594088E-005  [kg/m2]
 27:  end storage:   1.310177109332794E-002
 27:  pre_h2oveg  1.314278585736063E-002
 27:  ENDRUN:
 27:  ERROR in FatesPlantHydraulicsMod.F90 at line 2808  

In both cases I'm using this E3SM branch with this FATES branch (with the change to FatesHydroWTFMod for the test case). I can try also testing with this PR and the hydro stability PR.

The region used in these runs is a lat/lon rectangle from 15 N - 25 S and 30 W - 85 W so it includes very dry regions outside the Amazon basin where FATES-Hydro is unhappy.

Unfortunately I also get this error with a domain restricted to the Amazon Basin and when including @xuchongang's hydro stability branch. Also, I forgot to mention that I am using fates_hydro_solver = 2 .

@glemieux
Copy link
Contributor

glemieux commented Mar 7, 2024

I should note that the regression tests for this were all b4b otherwise.

derecho location: /glade/u/home/glemieux/scratch/ctsm-tests/tests_pr1164-fates-correctbase

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

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

@glemieux glemieux added the stuck Pull request is currently waiting on resolution of an issue that came up during review or testing label Mar 19, 2024
@rgknox rgknox added PR status: Not Ready The author is signaling that this PR is a work in progress and not ready for integration. and removed stuck Pull request is currently waiting on resolution of an issue that came up during review or testing labels Mar 25, 2024
rgknox and others added 7 commits March 25, 2024 13:05
…defined variables, reduce minimum suction of cch and fixed an error check
Bug fix. Adrianna Foster identified that trunks should be removed from the patch-level variable sum_fuel before being used to calculate rate of spread and fuel consumption.
…luxes to soil layer fluxes when no conductance (trivial bug fix).
@rgknox rgknox removed the PR status: Not Ready The author is signaling that this PR is a work in progress and not ready for integration. label Apr 17, 2024
@rgknox
Copy link
Contributor Author

rgknox commented Apr 17, 2024

These changes seem to be adding some stability. I ran a 25 year smoke test on the F45 grid with ELM-FATES-Hydro. Vegetation was initialized as a cold start using the specified stem diameter (which I specified as the diameter at 75% maximum height) and closed canopy assumption, along with fixed biogeography and no-competition assumptions. All parameters were defaults, with the one exception being the use of the Picard hydraulics solver.

The model completed the smoke test. I show some mean fluxes and states from the final year of simulation. I interpret from the output that the hydraulics model continued to generate fluxes within the realm of realism through the final year of simulation. Note that contrary to the Picard solver, the 1D Taylor solve (default) did crash after several years, having trouble generating a solution within the error tolerance.
f45-hydrotest-biogeonocomp-yr23.pdf

@rgknox
Copy link
Contributor Author

rgknox commented Apr 17, 2024

I think this PR is ready for another round of review, and has changed significantly since the first round. I plan no further changes beyond what comes out of review.

@jennykowalcz
Copy link

Wow, this is great @rgknox ! When perlmutter comes back up I'll try rerunning my various failed hydro cases with this branch and hopefully close several of my hydro Issues.
Should we make the Picard 2D solver the default?

@rgknox
Copy link
Contributor Author

rgknox commented Apr 17, 2024

I'm in favor of making Picard the default. I'm curious if others have found similar things. I'll create a discussion!

@jennykowalcz
Copy link

Wow, this is great @rgknox ! When perlmutter comes back up I'll try rerunning my various failed hydro cases with this branch and hopefully close several of my hydro Issues. Should we make the Picard 2D solver the default?

This PR with fates_hydro_solver=2 allows the cases to run successfully. Like you found, with fates_hydro_solver=1, we still get "Could not find a stable solution for hydro 1D solve" errors in dry gridcells

@bchristo
Copy link

bchristo commented Apr 22, 2024 via email

@@ -44,7 +50,8 @@ module FatesHydroWTFMod
! elastic-caviation region


real(r8), parameter :: min_theta_cch = 0.01_r8 ! Minimum theta (matches ctsm)
!real(r8), parameter :: min_theta_cch = 0.05_r8 ! Minimum theta (matches ctsm)
real(r8), parameter :: min_psi_cch = -20.0_r8 ! Minimum theta (matches ctsm)

Choose a reason for hiding this comment

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

update comment to say 'Minimum psi'

!this%th_min = min_theta_cch
!this%psi_min = this%psi_from_th(min_theta_cch+tiny(this%th_max))
this%psi_min = min_psi_cch
this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max))

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.

Comment on lines +1358 to +1367
real(r8) :: pc
real(r8) :: kr
real(r8) :: dkr_dP
real(r8) :: sat_res
real(r8) :: alpha
real(r8) :: lambda
real(r8) :: Se
real(r8) :: deltaPc
real(r8) :: dSe_dpc
real(r8) :: dkr_dSe

Choose a reason for hiding this comment

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

Just a minor comment. My memory is a bit fuzzy on my Brooks-Corey model without looking it up; can comments on these variables be made to ease interpretability? And an overall comment citing a companion source for implementing the three main different regimes (full, smoothing, saturated).

Copy link
Contributor

Choose a reason for hiding this comment

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

I agree that it would be great to add comments on the units and meaning of each variables

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 created an issue here: #1193
I would just add the changes in here, but I'm not familiar with that model.

Copy link
Contributor

@xuchongang xuchongang left a comment

Choose a reason for hiding this comment

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

I finished my review. It looks great to me (on minimum water content and forced zero conductance) and I think it should resolve the issue of very dry conditions.

@@ -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

!this%th_min = min_theta_cch
!this%psi_min = this%psi_from_th(min_theta_cch+tiny(this%th_max))
this%psi_min = min_psi_cch
this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max))
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.

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.

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?

@@ -44,6 +50,9 @@ module FatesHydroWTFMod
! elastic-caviation region


!real(r8), parameter :: min_theta_cch = 0.05_r8 ! Minimum theta (matches ctsm)
real(r8), parameter :: min_psi_cch = -15.0_r8 ! Minimum theta (matches ctsm)
Copy link
Contributor

Choose a reason for hiding this comment

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

The comment should be minimum water potential (MPa)


! If the difference between psi and psi-min is greater than 10MPa
! just assume there is no effect of the minimum function (ie weight 0)
min_ftc_weight = exp(min_ftc_scalar*max(psi_min-psi,-10._r8))
Copy link
Contributor

Choose a reason for hiding this comment

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

is min_ftc_scalar this a parameter value read from parameter file that can be user-defined? can we put 10 MPa as a parameter for future testing?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

When psi gets to about 2 MPA from psi_min, the weighting factor gets incredibly small, the 10 MPa here is just to avoid numerical problems with that exponent function. We could test the min_ftc_scalar though.. I don't expect the choice of min_ftc_scalar to have any noticable impact on model output, but I could try a few different values to see if anything happens. We have the utility parameter fates_dev_arbitrary that could aid in this.

! low suction
call get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi)

ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight
Copy link
Contributor

Choose a reason for hiding this comment

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

I really like this implementation and I think it should resolve a good chuck of issues of overshooting

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

Comment on lines +1358 to +1367
real(r8) :: pc
real(r8) :: kr
real(r8) :: dkr_dP
real(r8) :: sat_res
real(r8) :: alpha
real(r8) :: lambda
real(r8) :: Se
real(r8) :: deltaPc
real(r8) :: dSe_dpc
real(r8) :: dkr_dSe
Copy link
Contributor

Choose a reason for hiding this comment

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

I agree that it would be great to add comments on the units and meaning of each variables


return(om_watsat,om_sucsat,om_bsw)

def CCHParmsCosby84T5(zsoi,om_frac,sand_frac,clay_frac):
Copy link
Contributor

Choose a reason for hiding this comment

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

put some citation would be great

Copy link
Contributor Author

Choose a reason for hiding this comment

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

added

@@ -2767,6 +2784,20 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
end do

! Second pass, apply normalized weighting factors for fluxes

Copy link
Contributor

Choose a reason for hiding this comment

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

Good catch for low water content. If the sum of weight is close to zero, then the layer-specific water uptake is proportion of the root x water content here. It might be better to explain more on the new disaggregation

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 updated the comment a little, should be more accurate. Does this look good @xuchongang ?

…ro, to enable first time-step issues in ctsm
@rgknox
Copy link
Contributor Author

rgknox commented May 1, 2024

tests ok, tests_0501-150922de
note: I also tested without the protections added in the photosynthesis routine for small leaf areas, when I disabled that there was b4b on all tests.
note: the two hydro tests in the fates suite now pass (previously fail)

@rgknox rgknox merged commit e262198 into NGEET:main May 1, 2024
1 check was pending
glemieux added a commit to samsrabin/CTSM that referenced this pull request Jun 4, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
6 participants