Huite
commented
2 months ago This is a cute idea.
The transformation can also be done ad hoc, wherever it's relevant, which has the benefit of things looking slightly less exotic. I think the main difficulty lies with the DiffEq.jl convergence machinery? You'd have to override those methods somehow.
This issue came out of the thought: what if the state is not the storage itself, but some value that can be transformed into the storage?
Say a state is $u$, and we compute the storage from that as $S = f(u)$ with some invertible function $f$. Then
$$ \frac{\text{d} u}{\text{d}t} = \left(\frac{\text{d}f}{\text{d} S}\right)^{-1}\frac{\text{d} S}{\text{d} t}, $$
where $\frac{\text{d} S}{\text{d} t}$ is the current water balance function.
Say we let $f(u) = \ln(e^u+1)$, which is the softplus function. This function is extremely close to the identity function for $u \gg 0$, but goes asymptotically to $0$ for $u \rightarrow -\infty$. Also
$$ \left(\frac{\text{d}f}{\text{d}S}\right)^{-1} = 1 + e^{-u}. $$
What effectively happens here is that for large storage there is hardly any change, but small storage values get spread out over all negative values of $u$, and a completely empty basin corresponds to $u = -\infty$. This has the nice property that we do not need to check whether $u$ is out of domain (range) because any value of $u$ corresponds to a positive storage.
For the implementation this would mean that the initial 'storage' state is $u_0 = f^{-1}(S_0)$, and in the water balance function we need to compute $S = f(u)$. This adds quite some (conceptual) complexity, but might behave quite nicely around basins that are drying out.