_set_idomain function; assign idomain = 0

[RdWitte]

The current version of _set_idomain does not (correctly) set the idomain to 0 for the top layer when the top layer has a thickness of 0.
For example; cells with a river level will not be calculated correctly if the idomain value is set to -1 at a node with the top layer having a thickness of 0.

Currently, i have added the following two lines in my model to account for this.
mask = np.logical_and(ds["thickness"][0] > -0.01, ds["thickness"][0] < 0.01)
ds["idomain"][0].loc[mask] = 0

Is this something that could be added in a future update?

[dbrakenhoff]

@RdWitte, thanks for posting this issue, I vaguely recall running into this issue as well.

If you feel comfortable, feel free to submit a Pull Request with those two lines added to the _set_idomain() function and we'll take a look and update the code. Otherwise, let me know, and we'll add it to our todo list.

[bdestombe]

Isn't it better to compare resistances instead of thicknesses, and ideally compare it to the resistances/conductances of the neighbouring cells? This arbitrary 1cm can have severe impact when given a low conductance.

[rubencalje]

Do you mean the method nlmod.layers.set_idomain?

This methods set ds['idomain'] to 1 (active) where the cell thickness is larger than 0, and -1 (vertical pass through) where it is not. For the top and bottom cells this value should ideally be set to 0 (inactive), but I did not encounter any problems by setting it to -1. Can you explain a bit more what is causing the error?

However, it it is better to set the values of the upper and lower inactive cells to 0. set_idomain is called frequently, so it should be some fast numpy logic.

I do not think your lines solve the problem, as it sets idomain of all cells with an absolute thickness smaller than 0.01 in the first layer to 0. The inactive cells could also be present in the second layer though, and I agree with @bdestombe that a value of 0.01 is quitrer arbitrary, I would just stick with a positive thickness (>0), like we do now.

[RdWitte]

@bdestombe I don't think the .set_idomain function has anything to do with resistances. The function sets nodes to inactive, active or non-existent.

I agree the value is quite arbitrary. I have done this because i found that Python sometimes calculates thickness values of X^10-30 instead of 0 and I'm using a minimum increment of 0,1 for the thickness of the layers.

@rubencalje In my model, some rivers heads were ignored (points = modelgrid, shape = area of surface water in model, colors = head)

In hindsight, the problem might've occurred because of the fact that some cells where given a very small thickness, and therefore given an idomain value of 1. Which led me to adjusting the idomain value for the top layer for these instances to 0.

[rubencalje]

This solved by PR #250.