Hillslope hydrology

[swensosc]

Description of changes

Add hillslope hydrology parameterization

Specific notes

Changes include multiple soil columns per vegetated landunit, additional meteorological downscaling, new subsurface lateral flow equations, and a hillslope routing parameterization.

Described in: Swenson, S. C., Clark, M., Fan, Y., Lawrence, D. M., & Perket, J. (2019). Representing intra-hillslope lateral subsurface flow in the community land model. Journal of Advances in Modeling Earth Systems, 11, 4044–4065. https://doi.org/10.1029/ 2019MS001833

Contributors other than yourself, if any: None

CTSM Issues Fixed (include github issue #): None

Are answers expected to change (and if so in what way)?

• FSDS answers change due to rounding differences, since the history field now uses a column-level variable instead of a gridcell-level one. (Confirmed using this branch.)
• The oldhyd test changes answers due to the removal of the origflag parameter, which was intended to replicate CLM4 behavior.

Any User Interface Changes (namelist or namelist defaults changes)? Adds namelist parameters:

• use_hillslope: Toggle to turn on the hillslope model
• downscale_hillslope_meteorology: Toggle to turn on meteorological downscaling in hillslope model
• use_hillslope_routing: Toggle to turn on surface water routing in the hillslope hydrology model
• hillslope_head_gradient_method: Method for calculating hillslope saturated head gradient
• hillslope_transmissivity_method: Method for calculating transmissivity of hillslope columns
• hillslope_pft_distribution_method: Method for distributing PFTs across hillslope columns
• hillslope_soil_profile_method: Method for distributing soil thickness across hillslope columns

Testing performed, if any:

• Twenty year simulations to check for water and energy balance errors
• aux_clm on Derecho and Izumi is bit-for-bit identical to ctsm5.1.dev169 aside from the changes mentioned above.

[billsacks]

@swensosc - there were a few minor cleanup things I noticed while doing this review that I just went ahead and did myself because it was easier than commenting on them.

• [x] See https://github.com/swensosc/ctsm/pull/5 and please merge if you're satisfied with the changes there.

[billsacks]

I ran testing on this branch with comparisons against #1738 (with #1738 rebased onto dev090 and the cdeps external change backed out on the hillslope branch). Most tests passed, and all tests that passed are bit-for-bit with the changes in #1738 .

Failures that we'll need to resolve at some point are:

• MCT tests (I have an inline comment about this)
• Fortran unit tests
• Field lists differ, but that's probably expected

[swensosc]

I have it in initCold, but that must not have worked (?). I'll take those lines out, and try again.

On Thu, May 26, 2022 at 8:25 AM Bill Sacks @.***> wrote:

@.**** commented on this pull request.

In src/biogeophys/WaterStateType.F90 https://github.com/ESCOMP/CTSM/pull/1715#discussion_r882726973:

• if (flag == 'read' .and. .not. is_restart()) then

•   this%stream_water_lun(bounds%begl:bounds%endl) = 0._r8
  

• end if

initCold should handle that situation (initCold is done in all cases, even if there is a restart file). We occasionally need some extra backwards compatibility code, but in those situations you should also be checking .not. readvar – note that here you are resetting to 0 even if the variable was on the restart file. But in this case, I'm thinking that initializing it to 0 in initCold should be enough.

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#discussion_r882726973, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57D5CEGVXIPBAHOEIKLVL6CU3ANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[swensosc]

I don't understand; in the second if block, where will

this%topo_col(c) be set if the first if block isn't activated?

On Fri, Jul 15, 2022 at 10:49 AM Bill Sacks @.***> wrote:

@.**** commented on this pull request.

In src/main/TopoMod.F90 https://github.com/ESCOMP/CTSM/pull/1715#discussion_r922342392:

•      if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then
  

•         this%topo_col(c) = atm_topo(g) &
  

•              + (col%hill_elev(c) - mean_hillslope_elevation(l))
  

•         this%needs_downscaling_col(c) = .true.
  

•      else
  

•         this%topo_col(c) = atm_topo(g)
  

•      endif
  

  The change you made at the end of UpdateTopo is getting closer, but I think it still isn't quite right (sorry for my vague explanation before): I think you need to move the hillslope block of code here (i.e., the conditional block starting with if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then) to outside of the this%needs_downscaling_col(c) conditional – so that the hillslope adjustments are made to all hillslope columns, whether or not needs_downscaling_col is true or false at this point.

  Also, the comment right above this loop is wrong now: it says that it it shouldn't matter that we set topo_col = atm_topo, but now it DOES matter, because it's needed for the later hillslope adjustment. Can you please fix that comment?

I suggest the following code to address these two things:

! This could operate over a filter like 'allc' in order to just operate over active

! points, but I'm not sure that would speed things up much, and would require passing

! in this additional filter.

do c = bounds%begc, bounds%endc

 if (.not. this%needs_downscaling_col(c)) then

    ! For any point that isn't already set to be downscaled, set its topo value to the

    ! atmosphere's topographic height. This is important for the hillslope block

    ! below. For non-hillslope columns, this shouldn't matter, but is useful if

    ! topo_col is written to the history file.

    g = col%gridcell(c)

    this%topo_col(c) = atm_topo(g)

 end if



 if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then

    l = col%landunit(c)

    this%topo_col(c) =  this%topo_col(c) &

         + (col%hill_elev(c) - mean_hillslope_elevation(l))

    this%needs_downscaling_col(c) = .true.

 endif

end do

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#discussion_r922342392, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57AGS67P4APKAEY5XDLVUGJDHANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[billsacks]

I don't understand; in the second if block, where will this%topo_col(c) be set if the first if block isn't activated?

Ah, yes, that's a key point that I probably didn't explain: If needs_downscaling_col was already set coming into here, then that implies that topo_col would also already have been set by update_glc2lnd_topo. The key here is that we want to use that topo_col as a starting point in this case rather than using the atmosphere's topo value as a starting point.

[wwieder]

@billsacks I know you're juggling too many things, but I wonder where we are with this PR?
@swensosc What other code or support is needed to get a version of HH available for use in CTSM?

[billsacks]

Sorry, this keeps being near the top of my list but not quite at the top. Realistically, given my juggling of ESMF and various CTSM things, I imagine being able to take this on in a month or so. Back in July, the main remaining things to do were to wait for @swensosc to address my final review points and then introduce a new system test covering hillslope. @swensosc can you confirm that you have addressed these review points as far as you know?

[wwieder]

Hello, I wanted to circle back to the game plan for this PR, which has been open for over a year now. It seems like @billsacks has other activities that take priority over this work (no judgement). From the limited applications that @swensosc has helped walk me through the HH model enables science that is ready to go. Is there harm in bring this to main, realizing that like NEON, MIMICS and other new capabilities we'll continue making improvements in the code as we exercise its capabilities more?

[billsacks]

@samsrabin and @swensosc - Recording a few things we talked about this morning:

Some to do's:

• [x] We should add some error checks to prevent running hillslope with some incompatible options. The options I can think of that are incompatible with having multiple vegetated columns on the landunit are:

  • use_init_interp = .true.
  • Any transient landcover change: see https://github.com/ESCOMP/CTSM/pull/1249#issuecomment-781554388 for details

• [x] We should add at least one system test that exercises hillslope hydrology. @swensosc will provide a surface dataset for this - probably an idealized dataset to avoid issues that can arise with realistic datasets. We should at least have one global test. @swensosc should provide guidance regarding a minimum reasonable test length to exercise most / all of the hillslope code. If we need multiple years to trigger some aspects of hillslope behavior, then we should consider doing a short global test together with a longer single-point test.

• [ ] We hope to be able to demonstrate that this PR doesn't change answers. However, answer changes are expected in the system test that uses the oldhyd test due to the removal of origflag. To verify that the changes in this PR don't change answers for other aspects of oldhyd, we should do a test from master with origflag removed from the oldhyd testmod and verify that this branch is bit-for-bit with that one-off.

Some additional notes:

• One possible answer-changing modification is a change of sat_lev from being set at runtime to set as a parameter, as 0.9 (in SoilHydrologyMod.F90). I could imagine roundoff-level changes if there is a difference in the precision treatment of this setting now. So if there are unexpected answer changes, a first thing to try would be to back out that change.
• Some notes for the ChangeLog, in addition to documenting the main hillslope hydrology changes in this PR:
  • Fixed logic for writing 3d time constant history variables to first history tape so that they indeed only write to the first history tape.
  • Removed support for origflag == 1
  • oldhyd test will change answers due to removal of origflag = 1

[mvdebolskiy]

Just a note. Since excess ice PR has been accepted:

• all the columns will get the same amount of excess ice within a grid cell when initialized with a stream.
• If there are changes in energy balance within the soil column caused by lateral water fluxes column elevation will change during runtime due to different subsidence.

[samsrabin]

Thanks for the heads up, @mvdebolskiy! It seems to me—keep in mind I just started thinking about hillslopes yesterday—that the first behavior is correct. The second one seems physically realistic and so "correct," but is it desirable? And/or, will we need to rework some code to handle it?

[swensosc]

Both of these are fine. The first could be changed if one had a need and the second is desired; we would like to see different evolution in response to differences in e.g. north/south aspect.

On Fri, Aug 18, 2023 at 7:47 AM Sam Rabin @.***> wrote:

Thanks for the heads up, @mvdebolskiy https://github.com/mvdebolskiy! It seems to me—keep in mind I just started thinking about hillslopes yesterday—that the first behavior is correct. The second one seems physically realistic and so "correct," but is it desirable? And/or, will we need to rework some code to handle it?

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1683947539, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57GIHM4KDGY5AXDZR3DXV5W7XANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[mvdebolskiy]

Though maximum amount of subsidence, considering current initial volume is less than 1 m. So, it might be ok to not include elevation update for now.

[samsrabin]

I've completed aux_clm testing with the FSDS output change reverted and most everything is OK. Exceptions:

• I couldn't run FSURDATMODIFYCTSM_D_Mmpi-serial_Ld1.5x5_amazon.I2000Clm50SpRs.cheyenne_intel because the parent tag, ctsm5.1.dev130, doesn't include a fix that was needed to get that test to work for me.
• SMS_D_Ld1_P48x1.f10_f10_mg37.I2000Clm45BgcCrop.izumi_nag.clm-oldhyd had true diffs in lots of variables. Note the oldhyd specification; I'm not sure if this is something we really care about answer diffs in. The normalized RMS errors are pretty large, though.
• SMS_D_Ld1_Mmpi-serial_Vmct.f45_f45_mg37.I2000Clm50SpRs.izumi_gnu.clm-ptsRLA fails to build.

This last one is weird. A nearly equivalent test, SMS_D_Ld1_Mmpi-serial.f45_f45_mg37.I2000Clm50SpRs.izumi_gnu.clm-ptsRLA, built fine. The only difference is that the failing one specifies _Vmct. All tests with that specification but the Intel or nag compilers build fine; it's just the combination of _Vmct and the gnu compiler that has the problem.

[swensosc]

For your second, point does oldhyd activate origflag? If so, I removed origflag, so one would expect diffs there.

On Fri, Sep 22, 2023 at 9:30 AM Sam Rabin @.***> wrote:

I've completed aux_clm testing with the FSDS output change reverted https://github.com/samsrabin/CTSM/commit/def9c20c54a3ec9a51deae0d8668764890579e39 and most everything is OK. Exceptions:

• I couldn't run FSURDATMODIFYCTSM_D_Mmpi-serial_Ld1.5x5_amazon.I2000Clm50SpRs.cheyenne_intel because the parent tag, ctsm5.1.dev130, doesn't include a fix that was needed to get that test to work for me.
• SMS_D_Ld1_P48x1.f10_f10_mg37.I2000Clm45BgcCrop.izumi_nag.clm-oldhyd had true diffs in lots of variables. Note the oldhyd specification; I'm not sure if this is something we really care about answer diffs in. The normalized RMS errors are pretty large, though.

SMS_D_Ld1_Mmpi-serial_Vmct.f45_f45_mg37.I2000Clm50SpRs.izumi_gnu.clm-ptsRLA fails to build.

This last one is weird. A nearly equivalent test, SMS_D_Ld1_Mmpi-serial.f45_f45_mg37.I2000Clm50SpRs.izumi_gnu.clm-ptsRLA, built fine. The only difference is that the failing one specifies _Vmct. All tests with that specification but the Intel or nag compilers build fine; it's just the combination of _Vmct and the gnu compiler that has the problem.

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1731617035, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57EDLPLAX3F5KRFIZOLX3WVHVANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[samsrabin]

Yup, that could well be the reason. Why exactly did you remove origflag, and could it be added back to this PR?

[swensosc]

It's cleanup that we wanted to do. So 'no' to adding it back.

On Fri, Sep 22, 2023 at 9:47 AM Sam Rabin @.***> wrote:

Yup, that could well be the reason. Why exactly did you remove origflag, and could it be added back to this PR?

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1731645625, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57GSVERV3IINDAW2U6LX3WXK3ANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[swensosc]

I think that means that the test for origflag could then be removed with this pr?

On Fri, Sep 22, 2023 at 9:47 AM Sam Rabin @.***> wrote:

Yup, that could well be the reason. Why exactly did you remove origflag, and could it be added back to this PR?

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1731645625, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57GSVERV3IINDAW2U6LX3WXK3ANCNFSM5TV7P5HQ . You are receiving this because you were mentioned.Message ID: @.***>

[samsrabin]

Removing the oldhyd test would also mean there's no test with snow_cover_fraction_method = 'NiuYang2007'; not sure if that's something we care about. The other namelist modifications are covered in no_subgrid_fluxes.

[samsrabin]

Unfortunately simply reverting the "remove origflag and fracice" commit (84f14ff) didn't resolve the diff.

[billsacks]

Thanks @samsrabin !

Regarding the oldhyd test: see my suggestion in https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1682969486 -- running a one-off baseline from current master where we remove the origflag setting in the oldhyd test; I believe the expectation is that this branch should be bfb with master if you remove the origflag setting on master. If that isn't bfb, it might need more investigation.

Regarding the mct test: I'd say let's drop it, since we are no longer supporting mct.

[samsrabin]

Ahh, I missed that suggestion. Will do!

[samsrabin]

It passes!

[samsrabin]

@billsacks, re:

We should add some error checks to prevent running hillslope with some incompatible options. The options I can think of that are incompatible with having multiple vegetated columns on the landunit are:

• use_init_interp = .true.
• Any transient landcover change: see add namelist option to enable multiple columns having a single pft #1249 (comment) for details

@swensosc mentioned that he thought use_init_interp might not be incompatible after all. It basically would give same value to all columns in a gridcell, which should at least work and is probably better than needing a full cold start.

What do you think?

[samsrabin]

@swensosc, thanks for addressing Bill's comments! I noticed you checked off a couple that don't seem to have code changes associated with them, though. Could you please clarify? Specifically, this one and this one.

[swensosc]

As bill suggested, I removed the call to HillslopeDominantPft, but I forgot to remove the routine itself. I just pushed that additional change.

On Fri, Oct 6, 2023 at 1:27 PM Sam Rabin @.***> wrote:

@swensosc https://github.com/swensosc, thanks for addressing Bill's comments! I noticed you checked off a couple that don't seem to have code changes associated with them, though. Could you please clarify? Specifically, this one https://github.com/ESCOMP/CTSM/pull/1715#discussion_r1316189897 and this one https://github.com/ESCOMP/CTSM/pull/1715#discussion_r1316191440.

— Reply to this email directly, view it on GitHub https://github.com/ESCOMP/CTSM/pull/1715#issuecomment-1751306521, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGRN57D3RNS6IGUR5VKN533X6BLRRAVCNFSM5TV7P5H2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCNZVGEZTANRVGIYQ . You are receiving this because you were mentioned.Message ID: @.***>

[billsacks]

@swensosc mentioned that he thought use_init_interp might not be incompatible after all. It basically would give same value to all columns in a gridcell, which should at least work and is probably better than needing a full cold start.

What do you think?

To me, the most critical requirement of init_interp is that, if you interpolate an initial conditions file to "itself" – i.e., a compatible configuration where you could use that file without running init_interp – then you should get identical answers to a run where you didn't use init_interp. This requirement is checked via the LII test. Currently, I'm pretty sure this would be violated with hillslope: all output columns for a given grid cell would come from the first (or maybe last) vegetated column in the grid cell.

This may sound totally theoretical, but I feel like this is important in practice: It seems like init_interp gets triggered more often than you might imagine, e.g., because it's left on accidentally or because of some change in memory allocation that requires running init_interp even though the resolution is the same. The latter happens when going from a spinup to transient case, for example... though I realize this is not yet an issue for hillslope, since hillslope doesn't support transient.

However, I understand Sean's rationale for wanting to enable this. Either way you need communication: Either (1) you make it an error to try to do init_interp with hillslope and give an informative error message telling users why this is disabled and what to do to allow it (possibly commenting out this error check? -- it feels messy, but is simple); or (2) you need to somehow communicate that, while this works, it doesn't work to the same level of robustness as we try to ensure for init_interp in general. My personal preference is for (1), but I'm happy to leave the decision up to you (@swensosc and @samsrabin ) and others in the CTSM group.

[wwieder]

@samsrabin if this is just a test or two away from being merged I suggest we move this PR ahead of the other crop, tillage, etc work that you're doing. Does this make sense to you?

[samsrabin]

@wwieder Sure, sounds good.

[samsrabin]

@swensosc When running my test in debug mode, I get the following error:

66:forrtl: severe (194): Run-Time Check Failure. The variable 'soilhydrologymod_mp_subsurfacelateralflow_$C_DST' is being used in '/glade/u/home/samrabin/ctsm_hillslope_hydrology/src/biogeophys/SoilHydrologyMod.F90(2320,37)' without being defined

That line is here; the issue seems to be not that it's undefined but rather that it's uninitialized.

Update 2024-01-12: This is resolved as of c7c6602.