Calculate tracer concentrations internally

[evetion]

Fixes #1105

Very much a work in progress. Currently calculates the concentration of hardcoded tracers, using no user input from tables (although that stage has been set). Ever since the rework from Bart in #1819 I felt calculating tracers would be simple to implement.

This uses the intermediate flows from update_cumulative_flows to update the mass balance in all basins using straightforward addition/subtraction for each solver timestep. As such it should be more accurate than Delwaq (which uses the mean flows at the saveat interval). I do see some notable differences from Delwaq which needs investigating. Continuity tracer looks good though, so I'm slightly optimistic.

My initial attempt was to add this calculation to the save_flow callback, at the same time interval we export to Delwaq, but since this only does addition over large timesteps (whereas Delwaq does integration) the results were quite off (especially in the beginning, hitting negative concentrations). The end results looked similar to this approach. I also looked into adding this stuff instead to the water_balance! methods, but I'm out of my depth in how the RHS works (in terms of sciml, caches, etc.).

• Generic parsing of user defined concentration input
• Decide on either using interpolation for temporal input or update (discretely) via callback: update: Use discrete callback for updating concentrations
• Hook up use_evaporation (and come up with a better name, evaporate_mass?)
• Compare to Delwaq and explain differences (test smaller saveat timesteps)
• Added missing UserDemandConcentration, enabling temperature or other discharges from such nodes.
• edit: In new issue UserDemand should return mass, not swallow it.

[evetion]

This seems good to test with toy models now, everything is hooked up and running. HWS model runs fine. In testing the new testmodel I found I had to calculate inflows first and recalculate the concentration before outflows, otherwise the inflow/outflow could be reversed and result in negative mass.

The new testmodel switches a linear resistance on and off based on the concentration in the basin 🎉. However, I do see some weird effects, like consistent negative flows in both the turned off Linearresistance and Terminal, and the fact that the concentration goes much higher than the control should allow for. I guess I made an unrealistic model.

PS Is it correct that a linearresistance can be switched on/off based based on control, but that setting the resistance very high, or flow to 0 does not work?

[visr]

PS Is it correct that a linearresistance can be switched on/off based based on control, but that setting the resistance very high, or flow to 0 does not work?

No that should also work.

[evetion]

So if you make a small toy model with huge flows and a control that acts like a pinball machine, results are bad (and at one time lead to a timeout (2h+) of the CI with the new tolerances). The current model can be improved for sure (probably by fixing the discrete control and introducing a levelboundary instead of the terminal, so the salt water is flushed better), but it now tests the different callbacks, and triggers the control correctly.

[evetion]

There seems to be a good trend agreement between Ribasim and Delwaq, and the saveat timestep does not matter too much (and choosing 8 hour would solve any differences). Ribasim underestimates/Delwaq overestimates some concentrations, which I'm looking into in #1901.

[SouthEndMusic]

In general it looks nice, I had just a bunch of code quality comments. I will look into the assumption the piecewise constant concentrations in time.

[SouthEndMusic]

You are indeed cutting some corners numerically by assuming that concentrations of flows are constant over a timestep. This could be improved, but then we might get the same fiddliness that we had when trying to improve flow integration before the state refactor. A structural improvement in terms of accuracy would be to add the substance masses to the states, but let's only go there if we need to go there.

[Huite]

You are indeed cutting some corners numerically by assuming that concentrations of flows are constant over a timestep. This could be improved, but then we might get the same fiddliness that we had when trying to improve flow integration before the state refactor. A structural improvement in terms of accuracy would be to add the substance masses to the states, but let's only go there if we need to go there.

What's common in most hydrological software is to use a full blown solver for integration, but to fully separate the solute transport equations from the flow equations, or in DiffEq terms, to create a secondary ODEProblem.
The primary reason to do this is that time scales may be radically different, and it's a shame to solve the nearly same flow problem over and over. As long as we do not consider density effects, the flow "physics" do not depend on solute concentrations directly -- although PID controllers maybe complicate things -- only through control measures?

[evetion]

There seems to be a good trend agreement between Ribasim and Delwaq, and the saveat timestep does not matter too much (and choosing 8 hour would solve any differences). Ribasim underestimates/Delwaq overestimates some concentrations, which I'm looking into in #1901.

Now that I have a fix for that issue, I assume that the underestimation is because cyclic UserDemands (from/to the same Basin) correctly return concentrations in Delwaq (the flows are merged, so only 1-return% is removed), whereas we don't here (issue #1894).

[visr]

Is it easy to check this by setting all UserDemands to 0? Would be nice to see the lines from #1849 (comment) on top of each other.

[SouthEndMusic]

For reference:

[SouthEndMusic]

As long as we do not consider density effects, the flow "physics" do not depend on solute concentrations directly

As long as we do not do PID control or continuous control based on substances, the flow and substance transport problems are not hopelessly intertwined. When discrete control happens based on concentration, we can not fully see it as 'offline coupling'. so we would have to alternate steps of the storage and the substance problems. Still then there are numerical simplifications, see my whiteboard post above.

Handwavingly, the solver knows more about what flows do in between timesteps than what we get in the output, so solving both problems simultaneously can make use of that.

[evetion]

Is it easy to check this by setting all UserDemands to 0? Would be nice to see the lines from #1849 (comment) on top of each other.

Great minds think alike 🎉

[evetion]

You are indeed cutting some corners numerically by assuming that concentrations of flows are constant over a timestep. This could be improved, but then we might get the same fiddliness that we had when trying to improve flow integration before the state refactor. A structural improvement in terms of accuracy would be to add the substance masses to the states, but let's only go there if we need to go there.

I'll happily admit to corner cutting 😎. However, and maybe @Huite can chime in, I feel that concentrations can't change that much over a timestep. They're very much related to flow itself, and the continuity (flow * 1) tracers mass balance is identical to the water balance.

As long as we do not do PID control or continuous control based on substances, the flow and substance transport problems are not hopelessly intertwined. When discrete control happens based on concentration, we can not fully see it as 'offline coupling'.

Apart from some UserDemand problems, I see very similar results to Delwaq now, so this approach would at least be fine for fraction plots (no control at all). After fixing UserDemand, we should start testing model stability with control. Maybe something for @Fati-Mon at a later stage?

[Huite]

@evetion

However, and maybe @Huite can chime in, I feel that concentrations can't change that much over a timestep. They're very much related to flow itself, and the continuity (flow * 1) tracers mass balance is identical to the water balance.

No, this is not the case. You can convince yourself with a simplified, somewhat extreme example: think of a steady flows (dQ/dt = 0 and dS/dt = 0) for over a long time (say a year); the flow solver has no reason to (adaptively) subdivide time. Assume your initial tracer state is 0, and has concentration 1 on the incoming flows. What happens?

@SouthEndMusic

so we would have to alternate steps of the storage and the substance problem

Indeed.

the solver knows more about what flows do in between timesteps than what we get in the output

This is also correct, although here at least we do have a choice, unlike the storage <--> reconstructed flow conundrum, since that introduces a water balance inconsistency in the output.

It "merely" becomes a discretization choice in which we assume piecewise constant flows over a saveat timestep when we consider solute transport. I assume Delwaq does exactly the same, since it only receives integrated flows.
From a consistency and correctness point of view, incorporating everything in a single coherent system of equations is indeed best. But the increase in accuracy is likely marginal compared to the very bad effects it might have on performance (having to formulate and solve nearly identical flows over and over): see the difference between the saveat times that Maarten tested.

[evetion]

With the latest commit, this has probably also solved #1894 (will check in that issue).

[evetion]

However, and maybe @Huite can chime in, I feel that concentrations can't change that much over a timestep. They're very much related to flow itself, and the continuity (flow * 1) tracers mass balance is identical to the water balance.

No, this is not the case. You can convince yourself with a simplified, somewhat extreme example: think of a steady flows (dQ/dt = 0 and dS/dt = 0) for over a long time (say a year); the flow solver has no reason to (adaptively) subdivide time. Assume your initial tracer state is 0, and has concentration 1 on the incoming flows. What happens?

Concentrations in basins only directly next to boundaries are probably correct, the rest of basins are completely incorrect (not propagated to at all). Do such scenarios happen in our models though? My intuition for our (test)models (and hence my question) was that there are enough timesteps/discretization to actually propagate such changes. Do we allow for flows that would completely refresh basins (multiple times)? Does the saveat timestep trigger an integrator callback?

I've updated the docs to state our assumptions, feel free to suggest changes.

[Huite]

Do we allow for flows that would completely refresh basins (multiple times)?

Yes, this doesn't have any numerical or mathematical problems for flow.

Martijn and I used to run MT3D(MS) models, which simulates advective-dispersive solute transport. You could either use a low-order (stable) implicit scheme (with terrible numerical dispersion) or one of several explicit schemes; these all used criteria for adaptive time-stepping while holding water flow constant.

One of those criteria computes the refresh time of a cell, and sets the time step to the solute simulation to the smallest one encountered. Given the small time steps we currently encounter in Ribasim flow, those refresh times probably aren't a problem in general, but it's not guaranteed. That limitation is worth documenting for now.

My suggestion is to upgrade this implementation to another, separate, ODEProblem in due time.