I'm simulating temporal data of soil moisture $s(t)$ at a daily time scale by generating a random function with a compound Poisson process.

poi = CompoundPoissonProcess[1/6, ExponentialDistribution[1/15]]
data = RandomFunction[poi, {0, 10000}]

Where the occurrences of rainfall events are modeled using a Poisson process (and therefore have exponentially distributed inter-arrival times with a mean of 6 days in this case) and the depth of each rainfall event is also sampled from an exponential distribution with mean 15 mm/event. The RandomFunction simulates the process for 10000 days.

Between rainfall events (drydown), soil moisture decays through evapotranspiration through a function $\eta*s(t)$ where $\eta$ is the rate of evapotraspiration loss per unit ground are and in this case is 5 mm/day.

It is important to mention that, since our simulation is at a daily timescale, rainfall events lose all their temporal structure. This means that, mathematically, those correspond to sudden "jumps" or dirac deltas.

All that is left is to solve the differential equation $s'(t)=\frac{1}{soil\ volume}[f(t)-\eta*s(t)]$ where f(t) is defined as the forcing function (in other words, when the rainfall events occur). To generate f(t), I do this:

tdcum = data["Path"];
td = Differences[tdcum];
times = tdcum[[2 ;;, 1]];
depths = td[[;; , 2]];
f[t_] := Total[depths*DiracDelta[t - times]];

Here, the times variable stores the times at which the rainfall events occur and the depths variable stores the depth of each rainfall event. So I can create a function with DiracDeltas that will form the forcing of soil moisture. Something similar to this:$f(t)=1.81954 *\delta[-9998.13 + t] + 9.06213 *\delta[-9997.96 + t]+...$

The issue:

At this point I feel enough background has been given to now talk about the real issue. After soilving the differential equation sol = NDSolve[{ss'[t] == 1/400 (f[t] - 5 ss[t]), ss[0] == 0.4}, ss, {t, 0, 10000}]; and plotting at a certain intervals, I get cases like this:

s[t_] = ss[t] /. sol;
Plot[{s[t], 1}, {t, 4500, 4590}, PlotRange -> All, 
 AxesLabel -> {"times (days)", "s(t)"}]

enter image description here

Of course, the volumetric soil moisture should always be a number between 0 and 1. In some case, it rains more than the soil can admit and, in real life, the additional water gets transported through runoff.

I'd like to use the event detection capabilities of NDSolve to basically limit soil moisture to 1. However, it seems like this method struggles in detecting crossings of jumps as shown here:

sol = NDSolve[{ss'[t] == 1/200 (f[t] - 5 ss[t]), ss[0] == 0.4, 
    WhenEvent[ss[t] == 1, Print[t]]}, ss, {t, 4500, 4590}];

from which I get for my particular case shown in the figure above: 4555.09 and 4580.8 (2 crossings instead of the apparent 4). Those two times correspond to the times of down crossing (useless for my case) from above 1 to below 1. What I need is the time of up crossings to prevent soil moisture from exceeding 1. As you can see, WhenEvent overlooks up crossings committed by dirac delta "jumps".

Does anyone have a solution for this with or without WhenEvent?

  • 1
    $\begingroup$ Greetings! Make the most of Mma.SE and take the tour now. You have helped us to help you writing an excellent question. You can always Edit if improvable, to show due diligence, give more context, or improve your minimal working examples of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$
    – rhermans
    Apr 27, 2017 at 20:13

1 Answer 1


How about a solution with more WhenEvents? We can use WhenEvent to simulate the rainfall process directly within NDSolve. Then we can easily keep an eye on ss[t] to make sure it doesn't exceed 1 (here I chop it with Min):

sol = NDSolve[{ss'[t] == -5/200 ss[t],
  ss[0] == 0.4, tnext[0] == RandomReal[ExponentialDistribution[1/6]], 
  WhenEvent[t == tnext[t], {
    ss[t] -> Min[1, ss[t] + RandomReal[ExponentialDistribution[1/15]]/200], 
    If[ss[t] == 1, Print[t]], 
    tnext[t] -> tnext[t] + RandomReal[ExponentialDistribution[1/6]]}]},
  {ss, tnext}, {t, 0, 10000}, DiscreteVariables -> tnext];

Plot[ss[t] /. sol, {t, 0, 10000}, PlotRange -> All,
  AxesLabel -> {"times (days)", "s(t)"}, PlotPoints -> 400]

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.