hstepsize and nstepsize

[Hydrogeologist]

I’m using TimML to generate capture zones, and I’m seeing that the capture-zone shape changes noticeably when I modify hstepmax and nstepmax.

What is the recommended way to choose optimal values for these two parameters?
Should hstepmax be treated like a numerical tolerance (e.g., based on domain size), and should nstepmax simply be set large enough to avoid premature termination?

Any guidance or rule-of-thumb values would be appreciated.

Thank you!

[dbrakenhoff]

Hi @Hydrogeologist,

If possible, could you maybe share some figures showing the differences you are seeing? You're right about nstepmax. I thought TimML would tell you when a trace reaches nstepmax by printing a message but turns out it only does so if silent=False in w.plotcapzone().

You can always inspect the traces yourself by setting return_traces=True and metadata=True. With both of these set to true each trace is returned as a dictionary. The message key contains information about the termination of the pathline. For capture zones I guess you want your messages to say 'reached tmax'. So guidelines for nstepmax would be, just increase it until you're no longer seeing the reached nstepmax iterations message.

w = timml.Well(...)
# ... solve model etc.

traces = w.plotcapzone(
    nt=20,  # number of pathlines
    tmax=20 * 365.25,  # maximum time, days
    color="C1",
    return_traces=True,
    metadata=True
)
[i["message"] for i in traces]  # show termination messages

As for hstepmax, I'm not entirely sure about all the details as to how its used, but it represents the maximum horizontal step a pathline can take, so I think your analysis is correct and it should be based on your domain size. Smaller is more accurate but requires more computation time. So maybe a bit of trial and error is needed to find a good value for your problem. I guess once you find something that looks realistic, you could check by halving hstepmax and see if the result stays more or less the same.

Hope that helps a little. And maybe @mbakker7 has some more insights.

[mbakker7]

hstepmax is the horizontal step size (in meters or feet or whatever length units you are using). It will be shortened when the pathline enters a different layer, a different inhomogeneity, crosses an element, reaches an element, reaches tmax, moves vertically more than vstepfrac (as fraction of vertical thickness), or reaches the window boundary.
The user probably knows approximately what the pathline will look like and how accurate they want the pathline to be, and should pick hstepmax accordingly (the pathline will be a bunch of straight segments with length hstepmax).

A pathline will also terminate when nstepmax is reached, which is rarely what a user wants. nstepmax is mostly there so that the traceline won't go on forever when the user specifies inappropriate values or does something that wasn't expected by the modelere. It would probably be a good idea to throw a warning when nstepmax is reached, so that the modeler is aware of that. You can submit an issue on GitHub to add such a warning.

This can all be better explained in the docs. Maybe a 'HowTo' for pathlines and capture zones would be useful. Feel free to submit a notebook with such a HowTo if you want! Thanks for a good question!