Effect from ground water table depth seems unstable

[halla]

Hi, I'm playing around with the groundwater effect and getting unexpected results. Here is a plot of the AquaCrop_OSPy_Notebook_1.ipynb notebook example case, first season, with reduced precipitation, and running it with different constant water table levels. Any idea what is going on here? Is it a feature of the model or maybe something with the implementation?

Code to reproduce with the AquaCrop_OSPy_Notebook_1.ipynb notebook example:

import numpy as np

depths = np.linspace(1.0,5.0)
df_result = pd.DataFrame()
weather_data2 = weather_data.copy()
weather_data2['Precipitation'] = weather_data2['Precipitation']/10 # too much rain for ground water effect in the original
for d in depths:
    gw_model = AquaCropModel(SimStartTime=f'{1979}/10/01',
                      SimEndTime=f'{1980}/04/30',
                      wdf=weather_data2,
                      Soil=sandy_loam,
                      Crop=wheat,
                      InitWC=InitWC,
                      Groundwater=GwClass(WaterTable='Y',
                                       dates=[f'{1979}/10/01'],
                                       values=[d])
                    )
    gw_model.initialize()
    gw_model.step(till_termination=True)
    gw_model.Outputs.Final['water_table'] = d

    df_result = pd.concat([df_result, gw_model.Outputs.Final])

plt.plot(df_result['water_table'], df_result['Yield (tonne/ha)'])
plt.xlabel('water_table')
plt.ylabel('yield')

Maybe related, If I pick a depth of 2.142857 , (low yield case), one soil layer seems 'empty' from the start.

[thomasdkelly]

Hi @halla, yeah thats odd thanks for all the detail will take a look today

[thomasdkelly]

Running your code inside Notebook 1 i get this error when the model tries to initalize. Did you ever get this?

[/usr/local/lib/python3.7/dist-packages/aquacrop/initialize.py](https://localhost:8080/#) in read_groundwater_table(ParamStruct, GwStruct, ClockStruct)
    398             # if only 1 watertable depth then set that value to be constant
    399             # accross whole simulation
--> 400             zGW = df.reindex(ClockStruct.TimeSpan, fill_value=df["Depth(mm)"].iloc[0],).drop(
    401                 "Date", axis=1
    402             )["Depth(mm)"]

TypeError: value should be a 'Timestamp' or 'NaT'. Got 'float64' instead.

This is most likely down to the pandas version so will try and fix asap, but wondering how you got around it?

[halla]

Right, there was some problems with the latest version with pandas and pathlib if I remember right... I seem to be able to get the original online notebook working by downgrading pandas and restarting the runtime.

!pip install pandas==1.2.5

[thomasdkelly]

Great thanks, the version problems will be solved in the latest release, there are quite a lot of new structural changes to come over the next few weeks/month.

I have found the main error: Line 1107 just needed to be changed to for ii in range(compi+1):. https://github.com/aquacropos/aquacrop/blob/78ed597683028aa43261921ee73a0746700b24b8/aquacrop/initialize.py#L1106-L1112

After I make the change I get a much more reasonable output. There is still a small bump so will find out where that comes from

[halla]

Thanks for looking into this. The modification seems to fix the first gap but another remains (or maybe this is the 'bump'?)

Different soil types seem to produce a similar pattern, but Clay is quite special..

[thomasdkelly]

Yeah thats not great is it.

Found another of the same error which helps a little more:

https://github.com/aquacropos/aquacrop/blob/78ed597683028aa43261921ee73a0746700b24b8/aquacrop/solution.py#L367-L369

But still getting those spikes so will keep digging

[thomasdkelly]

Last thing for now:

So I think it may all be down to rounding / floating point errors. In this section of code dth is often very close to zero (e.g., <1e-5) but as it is just above zero we enter the if statement.

https://github.com/aquacropos/aquacrop/blob/78ed597683028aa43261921ee73a0746700b24b8/aquacrop/solution.py#L1754-L1759

By changing line 1755 to

dth = np.round(NewCond.th_fc_Adj[compi],3) - np.round(NewCond.th[compi],3)

The yield difference went from 4 t/ha to 0.3 t/ha.

[halla]

Thats right, I also remember seeing some strange floating values around. There is also a function math.isclose() for floating point comparisons, introduced in Python 3.5 https://peps.python.org/pep-0485/ , maybe worth considering.

[pacs27]

I have been working on the bug but I can't find it.

I have seen that using the data that appears in notebook 1 (SandyLoam soil) from 2.67 the bug disappears again (in that type of soil). The problem appears on day 140 of the simulation.

Here are some charts that I have made of the crop growth in the two scenarios. I leave them in case they can help you to know more about the problem.

[thomasdkelly]

I have found the bug mostly disappeares once adding the rounding edits above, did you add these when producing these figures @pacs27 ?

[thomasdkelly]

Maybe down the road we can add some kind of debugging mode that tracks more variables so that we can compare two different models more easily.

[pacs27]

I have found the bug mostly disappeares once adding the rounding edits above, did you add these when producing these figures @pacs27 ?

yes, those graphs are applying your change.

[pacs27]

Maybe down the road we can add some kind of debugging mode that tracks more variables so that we can compare two different models more easily.

Yes, it can be a good solution. Have you thought about how to do it?

[thomasdkelly]

Since all the calc functions are compiled it might would have to be a simple(ish) solution like creating a dictionary called info or logs at the start of each function. Then add all extra variables to that dictionary and return it at the end. And this info/logs dictionaries just get stored each day. Could run into memory issues though.

[pacs27]

I may have found the bug. I haven't compared it with the value given by the desktop version and it still gives strange values on the Clay type soil.

The problem comes because in the function check_groundwater_table which did not take into account the last compartment. A few lines before it subtracts 1 from the comp variable (I suppose to avoid having to put minus one when calling the array item) but then it doesn't add one in the range.

These are my results:

Sandy loam:

Clay:

SiltLoam:

Line with the bug: https://github.com/aquacropos/aquacrop/blob/78ed597683028aa43261921ee73a0746700b24b8/aquacrop/solution.py#L368

Should be:

for ii in range(compi+1):

[thomasdkelly]

Yeah thats not great is it.

Found another of the same error which helps a little more:

https://github.com/aquacropos/aquacrop/blob/78ed597683028aa43261921ee73a0746700b24b8/aquacrop/solution.py#L367-L369

But still getting those spikes so will keep digging

Is this not the same line

[pacs27]

wow.. I don't know how I haven't seen it.

[halla]

Hi again, any updates on this? Seems at least some changes made it to 2.2, right? I found some new 'interesting' behavior when exploring the soil moisture dynamics, at least some related to bottom layer and its relation to the watertable. I'm running the model with static weathers (eg, zero rain, zero C, constant ET) to isolate the soil behavior from plant growth. I could do some proper examples when I have more time.

[halla]

The specific problem I had was when using linear watertable 1...1.5m, with 1.2m soil . The bottom layers remain 'stuck under water' even after the watertable falls below 1.2m. As I workaround I shrinked the soil column to 0.9m to fix that.

(There's another interesting feature in the picture that the 'initial water content' is reset to field capacity on the planting date, is this intentional?)