Peter.

Thanks for your reply.

A few follow-up questions and comments.

1. You say that for a particular element, the "mean" conditions of the nodes on that element are used such that if I have an element with one node possessing a large inflow, then the in transfer coefficient will be used. I do not understand this since I have assigned a value of zero to the in transfer coefficient to all elements in the model. A non-zero out transfer coefficient is used for all elements, on all layers in the neighborhood of the drain. So even if there would arise an in-flow, q_in = -phi_in * (h_R - h) = 0 * (h_R - h) = 0.

2. I've tried using a maximum flux constraint of Qmax = 0. This works much better although I still see flow from the boundary into the model domain at very small (10-4 m3/d) rates. I can live with that. The problem with this approach is that the model runs VERY slow, for example, 63 hours to simulate 5 days. This is unacceptable when I need to simulate at least 60 days. Can you suggest anything to speed this up?

3. Have I defined the drain properly? I have defined h_R = z, h_min = z, both time constant.

If the head in the aquifer, h, is greater than the head assigned to the boundary (h_R), there should be out flow, that is, flow from the model domain into the boundary, the normal drain operation.

For the reverse condition, equation 14-7 of the White Paper says q = 0 when h_R(t) = h_min(t) AND h(t) < h_min(t). This is what I expect a drain to do, that is, become inactive with Q_in = 0 when h < h_R (= h_min).

What happens in the normal case since I have h_R(t) = h_min(t) AND h(t) > h_min(t)? Is the flux calculated by: q = -phi_out * (h_R - h) or does the contraint condition apply such that q = q_min = -phi_out * (h_R - h_min)?

Dwaine