Added PTF for estimating air entry pressure

[PeterNelemans]

Issue addressed

Fixes #28

Explanation

I defined an extra PTF for estimating the air entry pressure of soil, according to the equation from the paper in the issue. It's all part of the setup_soilgrids, so the changes are relatively small. Only downside is that the PTF is only valid for certain soils (percentage clay between 5-60%, sand between 5-70%). For these cells NaN is returned, which Wflow should handle by defaulting to -10 cm.

Checklist

• Updated tests or added new tests
• Branch is up to date with main
• Tests & pre-commit hooks pass
• Updated documentation if needed

I foresee two potential issues, that I would appreciate your thoughts on.

1. In workflows.soilgrids.py, the variable hb is first computed at the resolution of the source data, and then reprojected to the model resolution using the average. In order for this to work, the soil porosity thetas has to be kept in memory in the original projection, and exists temporarily next to the reprojected thetas (which I now renamed to thetas_out in this function). Thus, there is a bit of a memory overhang.

   I did consider computing hb at the model projection (and thus removing the need to keep thetas in memory), by reprojecting the sand and clay percentages (and using the already reprojected porosity). However, I think this could lead to undesirable behaviour, as the relation between soil composition and the air entry pressure is non-linear. Reprojecting the porosity, and clay and sand percentages may result in a air entry pressure that is not representative of the model grid cell. I prefer to compute the air entry pressure in the original projection, and then reproject the air entry pressure by taking the average.

2. If the soil doesn't meet the requirements for the PTF to be valid (5% < clay/sand percentage 60/70%), a NaN is returned. However, when reprojecting the air entry pressure from the original projection to the model projection, these NaNs are ignored. Therefore, the resulting PTF of a model grid cell may not be representative of that area (e.g. when only a single cell in the original projection meets the requirements, its air entry pressure will be used for the entire model grid cell). We can accept this, as it may only occur in rare cases, or we can use the default value (-10 cm) for these cells, as Wflow uses -10 cm as a default anyway. Or, for these specific cells, we can use some values from this report. E.g., for very sandy soil (>90%), we use -4cm, for very silty soil (>90%) we use -70cm, and for very clayey soil (>90%) we use -800cm (just taking some values out of the report). However, the report also mentions some issues with the original equations, so I don't know.

[PeterNelemans]

This is also my first PR through a fork, an anyway the first time working on hydromt_wflow. So, would be happy to receive any general feedback on the way of working, and maybe if I forgot something (e.g. updating the docs).

Also note that I didn't update any of the tests or examples, as that might be a lot of effort for somethign so small.

[shartgring]

Hi @PeterNelemans! It might be worth to check out https://deltares.github.io/hydromt/stable/guides/core_dev/contributing.html There you can find how to solve the linting issues

Also, some of the tests are failing. This is likely caused by changes in setup_soilgrids that you made. In the tests, a new model is build and compared to already existing model files (staticmaps, toml). These now give a mismatch so the reference model used in the comparison needs to be replaced

[PeterNelemans]

Thanks @shartgring for the tips!!! (Also for your help yesterday already with finding the .naming file). I will have a look

[PeterNelemans]

I've updated some stuff so that the tests should pass, but one still fails. Not really sure how to address that...

[vers-w]

Or, for these specific cells, we can use some values from this report. E.g., for very sandy soil (>90%), we use -4cm, for very silty soil (>90%) we use -70cm, and for very clayey soil (>90%) we use -800cm (just taking some values out of the report). However, the report also mentions some issues with the original equations, so I don't know.

I think a better alternative is to use the air entry pressure value based on soil texture from Clapp and Hornberger in Table 2.

[PeterNelemans]

Thanks for the fast reply @vers-w! So, would you opt for taking all air entry pressure values from this table, or just the ones that fall outside the range of the PTF? For clarity, everything within the red line is in the range of the PTF.

[vers-w]

Thanks for the fast reply @vers-w! So, would you opt for taking all air entry pressure values from this table, or just the ones that fall outside the range of the PTF? For clarity, everything within the red line is in the range of the PTF.

I would use the values from the table for the ones that fall outside the PTF range. And probably nice to add the Clapp and Hornberger lookup table also as an extra option for the user (so a user can opt for the PTF or Clapp and Horberger).

[PeterNelemans]

I would use the values from the table for the ones that fall outside the PTF range. And probably nice to add the Clapp and Hornberger lookup table also as an extra option for the user (so a user can opt for the PTF or Clapp and Horberger).

Good idea! I suggest to let the look-up table be the default for all soils, and the user has the option to enable the PTF, which would then be used only for soils within the range.

[PeterNelemans]

@vers-w, sorry it's been a while; busy with another project. Anyway, I still want to merge this at one point, and something similar is also in progress for hydromt_sfincs, for determining Green-Ampt infiltration.

I was looking into using this paper you mentioned:

I think a better alternative is to use the air entry pressure value based on soil texture from Clapp and Hornberger in Table 2.

Just out of curiosity, why the preference for this source?

And on that paper: unfortunately, they table only has 11 out of 12 soil classifications :( For whatever reason, they do not have the soil properties for silt. We could just use the values for silty loam, or we could use the values defined in appendix A.4 of this report. It's more recent anyway, and fully complete. But you might have a good reason for prefering Clapp & Hornberger? Hence my earlier question.

[LuukBlom]

I rebased your branch to the current main to resolve the merge conflicts!

So, before making any more changes, please git pull the changes @PeterNelemans

[PeterNelemans]

Hey @LuukBlom! Sorry for the late response, and thanks a lot for rebasing the branch! Much appreciated :)

Sorry I haven't worked on this branch for so long; I was swamped with another project, but will have time next week. I'll keep you posted!

[vers-w]

I was looking into using this paper you mentioned:

I think a better alternative is to use the air entry pressure value based on soil texture from Clapp and Hornberger in Table 2.

Just out of curiosity, why the preference for this source?

As discussed last week, also with Jaap, this paper is peer reviewed and very well known in hydrology/hydrological modelling. For silt we could indeed just use the values for silt loam. Besides that the report is less known, while scanning through it last week together, it was also not clear on how many (diverse) soil samples this was based.

[shartgring]

Hi Peter, nice progress! Some tips to help you out (I did not check your code in too much detail to see if all my comments apply):

• Do you need to change the sbom and pixi.lock? As no new packages are imported. You can revert the changes by taking the sbom and pixi.lock from main
• Pre-commit checks are not passing. You can run pre-commit by activating the pixi environment inside your hydromt_wflow git repo pixi shell and then running something like pre-commit run --all-files. If this doesn't work, perhaps first pre-commit init. It is also possible to use git via that same terminal, which means that when you git commit .... it will automatically run the pre-commit hooks
• Add the changes to the changelog. Are we introducing breaking changes (i.e. models created with the updated hydromt_wflow will be different from the original hydromt_wflow despite having the same config file)?. Then it is good to highlight this. Also if we introduce breaking changes, this means that the test models need to be updated
• We spend some time towards v1 to improve naming the staticmaps, see also Update the default names in staticmaps #393. So hb can be improved to fit the other layers

[PeterNelemans]

Thanks @shartgring!

1. No, I did not intend to change anything .pixi related. On my fork, a GitHub bot automatically made a PR (for my fork of hydromt_wflow) to update a bunch of .pixi related stuff. I figured it knows what it's doing, but I do see that indeed there are still small changes. Will revert those.
2. Did that already, but some tests are still failling. It's because the staticmaps produced by the updated workflow differs from the staticmaps in the example folder
3. Yes, will do, see also answer above.
4. Will do!

Thanks a lot for thinking along :) Will reach out probably at one point for some feedback