Python boundary condition on Richards_Mechanics process

Hello everyone,
I want to model the reservoir water level fluctuation and rainfall in the Richards_Mechanics process. To model the reservoir water level fluctuation and variable rainfall intensity I choose to use python boundary. But it seems not working with Richards_Mechanics and giving error saying "No d.o.f. found for (node=6389, var=0, comp=0). That might be due to the use of mixed FEM ansatz functions, which is currently not supported by the implementation of Python BCs. That excludes, e.g., the HM process."

This pull request documentation also says “Only works if all physical fields have the same ansatz functions”.

My question is that is it still possible to use python boundary condition in Richrds_Mechanics? If not what could be the possible option for that?
My rainfall and reservori water level data looks like this.

Trial Inputfile:Groundwater.zip (188.3 KB)

Thank you,
Laxman Kafle