R^2 computation for different periods and different simulation results when adjusting tmin

[lauramenezess]

Hello,

For model validation using Model.plot() we should add a tmax= 'the end of the validation period' and do NOT use a tmin='the beggining of the validation period', right? I was doing it this way to have the R² for the validation period, but I have just realized the last value of the simulated head doesn't match the first value of the simulated head when validating (something to do with the constant d I imagine?)... Is there any other way of having the R² only for the validation period in PASTAS or should I use another package?

Thanks.

[dbrakenhoff]

Hi @lauramenezess,

The best way to compute an R² for a specific period is to pass the tmin/tmax to ml.stats.rsq(tmin=<tmin>, tmax=<tmax>). You can also do it manually by computing the residuals between tmin/tmax (res = ml.residuals(tmin=<tmin>, tmax=<tmax>)) and passing that to ps.stats.rsq(res=res).

I'm a bit confused why the first value of the simulated head is different to the last value when applying different periods in in Model.plot(). Maybe you can upload your model as zip file (use ml.to_file("model.pas") and zip that file), or you could upload some results pictures of your model, so we can take a look? Pastas models have a warmup period, but this should be sufficiently long that any effects of the spin-up should not affect the simulation in the calibration/validation periods.

Kind regards,

Davíd

[lauramenezess]

Ok, thanks for your answer regarding the R².

For my other questoin I am not sure if I made myself clear. As an example I have used the calibration period from 2013 to 2020, and the validation period from 2020 to 2024. If I use Model.plot() with a tmax='2024' I have this:

But if I use Model.plot() with a tmin='2020' and a tmax='2024' I have this:

We have very different results for the validation period, so I assume we should not add a tmin while using Model.plot()...
Not really a question, it's more a consideration...

[dbrakenhoff]

Thanks for posting the plots @lauramenezess. Based on these results, I think this is a case where the warmup is too short, or more likely, the model isn't able to simulate the heads without a very long response time. I think if you plot the step response for each of the stresses, for one of them (if you have more than one) (ml.plots.step_response()) you will get a graph that has very long time-scale on the x-axis (longer than 10 years anyway).

That means the model isn't reaching a steady-state within that default warmup period of 10 years, and in that case, when you modify the date for which you want to compute the model result, you will end up with a different answer. There is a warning added to the ml.fit_report() at the bottom when the response times are long compared to the length of the calibration period. Do you see that warning?

Having said that, I do wonder why we do not add a warning specifically related to the warmup when this happens. It's probably worth opening an Issue for that.

I can't be sure exactly, but looking at your head observations, there seems to be a downward trend in the data. If the stresses you have added to the model do not explain this trend, then Pastas might start doing some "statistical acrobatics" to try to fit the data as well as it can. In this case, it picks a very long response, which allows the model to simulate that downward trend, which isn't physically realistic, but does provide the best fit. That would be my guess anyway...

If it makes sense for your situation, you could try adding a linear trend to the model to detrend the period 2013-2020 using ps.LinearTrend()? That trend then signifies some other unknown influence on the head. Maybe that will help reduce the response time of the other response function(s)?

[lauramenezess]

Thanks a lot @dbrakenhoff, it was exactly that. I had a very long response time for the pumping stress. For models with a time response lower than the calibration period I don't have this difference in the results when using a tmin for validation.

Actually I have many issues with my pumping stress starting that it's a fictitious one. We've tried to recreate a pumping signal, not really minding about the values, with the periods when the farmers would likely pump (june to september). But it's not always the case, sometimes they start earlier of after that, so we are still trying to adapt this considering the balance between precipitation, ETP and the plants need of water... What we know is that they are pumping a lot more in the last 10 years than before, which may be causing this downward trend in the heads that is not completely explained by the balance between rain-ETP as we would expect... changes in temperature etc influence the plants need of water (and so the pumping periods) and we do not take this into account (yet).

I also wonder how could we differentiate between a good model considering all these stresses or simply a model that has many parameters to adjust so it will surely give a better result (right?)...
I've gone completely off topic, sorry!

Best regards,
Laura

[dbrakenhoff]

Hi @lauramenezess, I'm glad to hear that you've found the cause of changing simulations depending on tmin!

As for the problem you are looking into, it's an interesting and challenging one. We have in the past tried similar analyses to quantify the effect of unmeasured pumping for irrigation. For the time series you show above, it does not really have observations that directly point to pumping (i.e. no observation showing drawdown and recovery curves), so if that trend is caused by excessive pumping, it might be a diffuse and longer-term effect, which probably leads to Pastas estimating that long response time in that specific model. In this particular case fitting a linear trend might be the best you can do, but tying that downward trend to the lowering of groundwater levels to excessive pumping is obviously challenging. And not as strong as when including your fictitious pumping time series produces a reliable model.

As for determining whether your model is "good", that really depends on the goal. And you're right that including more parameters almost always improves the fit. In this article we proposed some criteria for checking our models, and maybe you can check some other Pastas publications on how they checked their models. A good first step is obviously a visual inspection and checking whether the results are plausible from a physical/hydrological point of view. Also splitting the observations in a calibration and validation period, like you did, is also a good idea to test if the model continues to perform well for data it was not calibrated on.

Anyway, good luck, and if you have any interesting results or more questions, feel free to share them here!

[lauramenezess]

Perfect!
I will soon share here some examples of Pastas application for the Alsatian aquifer, in eastern France.
Thanks a lot for your help!