Significant differences in model output since v1.0.4

[LuukBlom]

Found this issue when running the regression tests for hydromt wflow.

Given the same models, the v1.1rc1 and v1.0.4 releases produce very different results.

I ran the whole regression pipiline using hydromt_wflow's last stable release (v1.0.2) and with the current state of the release branch (v1.0.3.dev0)

see https://github.com/Deltares/hydromt_wflow/blob/main/tests/regression/README.md for an in-depth overview of the regression setup.

Reproduce

git clone https://github.com/Deltares/hydromt_wflow.git
cd hydromt_wflow
git checkout release/v1.0

$env:WFLOW_CLI="abs\path\to\v1.0.4\wflow_cli.exe"
pixi run regression-pipeline piave tests\regression\.runs-v1.0.4
pixi run regression-assert piave tests\regression\.runs-v1.0.4

$env:WFLOW_CLI="abs\path\to\v1.1rc1\wflow_cli.exe"
pixi run regression-pipeline piave tests\regression\.runs-v1.1rc1
pixi run regression-assert piave tests\regression\.runs-v1.1rc1

From the figures below, we can say that the hydromt_wflow changes are likely not causing the differences.
Note that the baseline metrics were generated with wflow.jl v1.1rc1 in both cases.

hydromt_wflow v1.0.2, wflow.jl v1.0.4

hydromt_wflow v1.0.3.dev0, wflow.jl v1.0.4

See Also

Deltares/hydromt_wflow#829 (comment)

[visr]

If I understood correctly, significant changes in results are expected from #730.
We still need to make some kind of plan for how to deal with that before we release.

[JoostBuitink]

If I understood correctly, significant changes in results are expected from #730. We still need to make some kind of plan for how to deal with that before we release.

This is indeed correct!

[JoostBuitink]

Adding to this issue to document our findings with the Meuse model.

The figure below gives a good summary of the differences between the same model schematisation and the difference in Wflow version. Ignore the 0.7 results, as this is a different model schematisation with different parameters than the 1.0 and 1.1 lines. In the calibration, we allowed for tuning of the soil_thickness, f, kv, fkh, rooting_depth parameters, with rather flexible adjustment factors, but none of the combinations where able to get close to the observations, specifically for the recession/baseflow periods. In general, the rising limb of the hydrograph looked rather good (maybe a bit too flashy), but of the discharge was dropping too quickly to too low values.

This is very consistent with the testing I had done when reviewing the PR (see below, ignore the first year due to cold state). Initially we could explain these results due the expected difference in subsurface flow, and rising/falling dynamics of the water table. We hypothesized that changing parameters as subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio would be able to compensate for this behaviour, but did not test this hypothesis.

However, after quite significant efforts in the calibration of the Meuse model (as described above), we learned that even with quite large parameter ranges, we were unable to get discharges to not drop rather quickly to zero.

[vers-w]

@JoostBuitink : could you share the complete Meuse model setup?

[JohnRushKucharski]

I wanted to add a note regarding what I think this issue requires:

1. An explaination of the difference. I think there is a hypothesis that this can be explained by changes in how field capacity is handled (I think this is documented in the 4th bullet under the unreleased changelog changed heading. However, it needs to be confirmed that this can indeed explain the differences in the result, and that these differences are intentional.
2. Documentation. Assuming the source of the change can be verified and is intentional the wflow documentation needs better document this change and its implication for users upgrading models. This document should whereever possible indicate the expected directionality, magnitude or pattern of changes in outputs and steps users may need to follow to recalibrate their models.

[atsiokanos]

Perhaps not so relevant for this issue as it concerns v1.1-v1.0.4 comps, but we recently had some discussions and did some analysis in Sava with @alimeshgi and @aweerts about which version to use. In our case, we had an old calibrated model that we wanted to continue using. And this is a relevant question as we have many projects with many models we want/(or suggest?) to update to the latest.

There were quite some differences between the versions which we were trying to explain. I think this is something that should be documented clearly, with explanations of the changes and their impacts whenever there are breaking changes in a release. Maybe we could also add some example hydrograph figs, e.g. for the Meuse.

Our conclusions:

• From v0.7.3 to v0.8.1 flows appear generally higher. Probably this has to do with the fact that we have less evaporation due to the difference in ET partitioning and how water is extracted from the soil.
• From v0.8.1 to v1.0.3 the water extraction is a bit more in line with how it was in 0.7.3 but the main and largest difference probably comes from changes in the infiltration capacity of the soil (parameter infiltcapsoil is not supported and infiltration is now determined from the actual saturated hydraulic conductivity at the soil surface using parameters kv_0 and kvfrac).

[vers-w]

There were quite some differences between the versions which we were trying to explain. I think this is something that should be documented clearly, with explanations of the changes and their impacts whenever there are breaking changes in a release.

Yes, I agree that it would be good to document the impact of process improvements (I think that is a better term than "breaking changes") that may require model re-calibration. We started documenting this for the next release in the changelog. Not sure if we should also add figures for a catchment to the docs.

[vers-w]

For the Meuse model shared by @JoostBuitink I did notice that the theta_fc parameter is provided as input:

What I did notice is that theta_fc is probably not based on soilgrids data, the pattern is quite different from theta_s (see figure below)
and theta_r, and the values are quite low (unrealistic). On what data is theta_fc based?

While I think it is fine to not tune theta_fc during calibration. I would say it is critical/preferred that theta_fc is consistent with the other theta parameters (based on the same input data, e.g. with soilgrids using a PTF).

[vers-w]

For comparison this is the default theta_fc parameter in Wflow:

[AquaPore]

While I think it is fine to not tune theta_fc during calibration. I would say it is critical/preferred that theta_fc is consistent with the other theta parameters (based on the same input data, e.g. with soilgrids using a PTF).

As a soil physicist, I would agree that introducing a new FC parameter into the Wflow model, which is not linked to the $\theta(h)$ function, would remove the beauty that Wflow is a semi-physical conceptual model. Instead, I propose allowing flexibility while remaining physical. For that, I propose introducing a new parameter, h_FC [L], the pressure at which FC is computed. h_FC is influenced by texture; water held at h_FC=33 kPa is mainly used to determine the θFC for clayey and loamy soils, whereas h_FC=10 kPa is mainly used for sandy soils (Libohova et al., 2018; Robertson et al., 2021). Thus, typically h_FC could vary between [10, 33] kPa.

With kindness,
Joseph A.P. Pollacco

[vers-w]

Some results below for about the same period (Figure in #1050 (comment)) for station Chooz. With the default theta_fc parameter and manual tuning of f (scale = 0.5) and soil_thickness (scale = 2.0) parameters the results are reasonable and discharge is not dropping to too low values. With setting theta_fc equal to theta_r (basically the v1.0 implementation except v1.1 makes use of fillable porosity) and using the same tuned f and soil_thickness parameters the baseflow is dropping to lower values, especially in the second year and peaks are generally higher.

The run with theta_fc == theta_r shows much lower net recharge values than the run with the default theta_fc parameter. Could be interesting to check this further (differences in evaporation saturated zone, capillary flux). Finally, could be useful to check results with v1.1 without the use of fillable porosity.

[vers-w]

Thanks for your input @AquaPore , I think what you describe could be a nice addition to hydromt_wflow, besides using a PTF this could be an additional option/transfer function. This way we keep the kernel and data pre-processing nicely separated. As a default in Wflow we use h_FC =10 kPa.

[AquaPore]

Thanks for your input @AquaPore , I think what you describe could be a nice addition to hydromt_wflow, besides using a PTF this could be an additional option/transfer function. This way we keep the kernel and data pre-processing nicely separated. As a default in Wflow we use h_FC =10 kPa.

Dear @vers-w ,

A great thanks, what I am proposing to remain physical is that theta_fc is no longer used as an input parameter and we remplace it with h_FC which can be adjusted.

With kindness,
Joseph A.P. Pollacco

[AquaPore]

Dear Wflow developers,

From a hydrological perspective, we agree that storativity computed from $\theta s - \theta fc$ is more physical. Nevertheless, it is shown that it does not always improve the predictions; therefore, I am questioning if the processes of Wflow can be improved and make Wflow more physical, by better taking advantage of the input soil hydraulic parameters:

Improvement 1

From a soil physics perspective, one of the limitations of Wflow is that θs does not decrease with depth, and therefore I suspect that the value of θFC is too high. I published several papers which relate (θs-θr) to the vertical saturated hydraulic conductivity Ksv, which I can share. Therefore, since subsurface flow is below the surface, then θfc which is related to θs will also need also to be reduced accordingly;

Improvement 2

We foresee that significant improvement would be made if we combine subsurface flow with groundatwer flow . A proposed solution would be where subsurface flow is computed for every cell except for the bottom cell (reserved for the groundwater module) and except for the top cell (reserved for overland flow).

It is important that subsurface flow would be computed from the horizontal unsaturated hydraulic conductivity instead of the horizontal saturated hydraulic conductivity.

With kindness,
Joseph A.P. Pollacco