Inconsistencies with Polder function

[AviadAvr]

Hi all, and thanks for the great effort put into this amazing project!

Although i've been getting pretty good results with using the Gamma function with modelling surface water as a stress, i've been curious to try out the Polder function, as I understand that it is meant to be designed specifically for that.

I've been getting quite inconsinstent results with the Polder function, especially when changing the length of my series (both the head observations and surface water levels).

a bit of background:

• I typically model short periods (1-12 months) of high frequency observation data (1 hour intervals), but daily meteo data.
• The well I analyzed in this thread was in a coastal region (Noordwijk, NL), 40m deep and the surface water level was sea water level (Scheveningen, NL). The distance between the well and the surface water measuring station was approx. 17km, but as the well is heavily affected by tidal influence, this seemed like a good pick (and the nearest station i could find).

At first I gathered data for 3 months (head and surface water) and 10 years of meteo data, then solved using either the gamma or polder functions for the surface water stress, while changing the duration of the solving period (lowering the tmax param) from 3 months ->2 months -> 1 month and then plotting the results on top of the head observations.

def makeSomePasta(head, meteo, SeaLvl, rfunc):

    if rfunc == 'polder':
        rfunc = ps.Polder()
    elif rfunc == 'gamma':
        rfunc = ps.Gamma()
    
    ml = ps.Model(head)
    rch = ps.RechargeModel(meteo.precip, meteo.evap, name='recharge', rfunc=ps.Gamma())
    ml.add_stressmodel(rch)
    
    surfacewaterlvl = ps.StressModel(SeaLvl, rfunc=rfunc, name='SurfaceWaterLvl', settings='waterlevel')
    ml.add_stressmodel(surfacewaterlvl)
    
    months_to_deduct = [0,1,2]
    
    months = []
    EVPs = []
    for i in months_to_deduct:
        tmax = str(datetime.datetime.now() - datetime.timedelta(days=30*i))[:10]
        
        ml.solve(freq='h', noise=False, tmax=tmax)
        
        EVP = ml.stats.evp()
        duration = 3-i
        
        months.append(duration)
        EVPs.append(EVP)
    
    return months, EVPs

Here are the results:

With the Gamma function, I get a consistent 92% EVP for the different durations, while with the Polder function I get ~10% EVP for the first 2 months, and then up to 94% for modelling all 3 months.

Furthermore, when I Increase the length of the entire series to 12 months and then change the solving period, I get completely different results for polder:

Again with the Gamma function we see pretty consistent results at around >90% for solving in 4 different durations.
Polder now shows low EVP for the first 2 months, and for modelling a duration of 3 months in total now only results in 35% EVP instead of the previous 94%, and peaks at 37% for the entire duration. I've also noticed Polder takes a much longer time to calculate.

Is this approach flawed? Or should I keep using Gamma for surface water levels?
I will be happy to provide the data series, if needed.

Thanks!

[dbrakenhoff]

Hi @AviadAvr,

First of all, glad you like using Pastas and thank you for sharing your experiences!

As for your question, that's a challenging one... As far as I know, the Polder function hasn't been as thoroughly tested or used as Gamma, so you might be running into some inherent challenges with that particular response function. One hunch I have is that the Polder function might be sensitive to the initial parameters. I'm not entirely sure how Pastas computes the initial guess, but given that changing the length time series leads to difference results despite the calibration period remaining the same, suggests that the initial guess is performed on the entire time series, and not the calibration period. The optimization might be getting stuck in some local minimum and seems quite sensitive to that initial guess? But this doesn't quite explain why the result in the 3-month run gets worse when the calibration period is shortened.

A couple of things you could investigate:

1. How do the initial parameters differ between the runs you've performed? Does setting them to the initial guesses of the one successful run improve the fit result in the other cases?
2. What fit/result do you get when you use the parameters from your one successful run with the Polder function in each of the other cases? Judging by the results, my guess is those optimal parameters should result in a decent fit in each of those cases... This should prove a better result is attainable.
3. Does using a different solver give you different results? E.g. ps.LmfitSolve().

Anyway, those are some thoughts I had when reading your post, maybe they provide some ideas on how to further investigate this result. If you have any new results or insights let us know here!