Small number protections to CCH water transfer functions

rgknox commented 9 months ago

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

[x] The in-code documentation has been updated with descriptive comments
[x] 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

Technical Note update:
User's Guide update:

Test Results:

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

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

FATES baseline hash-tag:

Test Output:

rgknox commented 9 months ago

This may address: #1154

jennykowalcz commented 9 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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

rgknox commented 7 months ago

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 commented 7 months ago

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 commented 7 months ago

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 commented 7 months ago

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

jennykowalcz commented 7 months ago

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 commented 7 months ago

Chiming in here. This is fantastic news and sounds to me like the Picard solver should be the default.

On Mon, Apr 22, 2024, 10:33 AM Jennifer Kowalczyk (Jenny) < @.***> wrote:

Wow, this is great @rgknox https://github.com/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

— Reply to this email directly, view it on GitHub https://github.com/NGEET/fates/pull/1164#issuecomment-2069913077, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADSXQZRRO2RVUUSOQBUKPTTY6UUUFAVCNFSM6AAAAABDVVGDSCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANRZHEYTGMBXG4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

rgknox commented 6 months ago

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)

NGEET / fates