Representing drains

I am having difficulty with representing drains/interception trenches in my model.

I am using a third-type boundary condition, h at the drain nodes is set to z, the node elevations, and a minimum constraint h_min set also to z, the node elevations. The maximum constraint is not set (+ INF). Assigning the BC and the contraint is done by joining to a polygon representing the location of the drain. This assignment is made on adjacent slices (slices 9 & 10 of an 11 slice model). I have assigned the OUT transfer coefficient to be very high (for an ideal drain) and the IN tranfer coefficient set to zero. These transfer coefficients are assigned to each element within the polygon on all layers. The elements within the drain polygon are assigned a Kx=Ky=Kz of very high (for an ideal drain) on all slices. The mesh is very refined in the region surrounding the drain.

My problem is, after the model is done and examine the results in the budget analyzer, I see that that during certain times (transient model) I have water moving from the drain into the domain approaching the same rate that water is leaving the domian into the drain. My expectation is that water should never flow into the model domain from the drains. When there is a preoblem with influx, the in-flowing nodes seem to be consistent, that is, the influx always seems to occur in the same locations.

The purpose of the model is to detrmine the sustained flow rate into the drain so this precludes using type-4 well boundary conditions.

What can I do to ensure that the drain never produces water?

Hi Dwaine,

The reason for the infiltration is probably the following: Transfer rates are set to elements, while the transfer boundary condition is set to nodes. For deciding whether to use the transfer rate in or out for an element, FEFLOW uses the mean conditions. So if at one of the nodes at the element there would be a little bit of inflow, and at the other node a high outflow would occur, the transfer rate out is used for that element. Unfortunately that also influences the condition at the inflowing node - leading to an inflow.  So the inflow usually occurs at the nodes where inflowing and outflowing conditions change.

The only way to avoid inflow completely is a flux constraint. However, that would interfer with your 'constrained by head' for the drain. So  there's no 'perfect' solution for that case, the procedure highly  depends on the model details and aims.

Cheers,
Peter

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

Dwaine,

I'm sorry for the delay, but we have been traveling for a while without much time.

1. The problem is the following: Imagine you have two neighboring nodes, and at one node you'd get a q_in, for the other one a q_out because of changing conditions between the nodes. Then for the elements between the two nodes only one transfer rate can be used (as this is an elemental parameter). If the average conditions are 'out', then for BOTH nodes the out transfer rate will be used, even if there's an inflow to the node.

2. Unfortunately, no. For applying the min/max flux constraints, in each time level a budget analysis has to be done for all the nodes with constrained conditions - and that takes time. Theoreticall there would be another option in case of a transfer condition by just looking at the gradient, but it hasn't been implemented yet (we have only used that in custom programming using the FEFLOW programming interface).

3. Yes, that's OK. And it should work well, not giving you infiltration. For the normal case (h(t) > h_min(t)) the calculation is q = -phi_out * (h_R - h).



Cheers,
Peter