3-d Fluid-flux BC

Hello,
I am working on transferring an existing 3D finite element regional flow model to FEFLOW (steady-state, flow only). I have successfully transferred my mesh (triangular prisms), material properties and my constant head boundary conditions. However, I am a little bit confused as to how to apply the recharge (Neumann BC).
In my initial model, the recharge, which is spatially variable, is defined per element face in mm/yr. On the other hand the Feflow Neumann BC is defined per node (L/T). The user manual specifies that to be valid, the BC “must be applied to all the nodes of the element face”.

My questions are:
What happens if two adjacent elements have different recharge values?
What is the best way to apply a spatially varying recharge?

My initial approach was to use arcmap to calculate equivalent recharge values per nodes using the intersection of the node’s thiessen polygons to the elements. I obtain a recharge value in mm/y which would apply to the area of the thiessen polygon.

However, once I applied my equivalent nodal recharge to my model and run the model. The model does not converge and the error is oscillating from one iteration to the other. The rate budget analysis gives recharge values that are 3 orders of magnitude larger than the value given in the mass balance of my original model. And also, I obtain both a considerable outflow and a considerable inflow through the Neumann boundaries.

I have double-checked the units and the sign of the applied parameters and they seemed to be consistent with the original model.

Thank you for the help!

Debora

Recharge is assigned using the material property "In/out flow on top/bottom" rather than a BC.

Hello,

Thank you Blair for the clarification. After applying the recharge as a material property, the rate budget does indeed resemble more the water balance from the previous model.

However, my model still fails to converge. (I set the convergence criterion at 1e-5 as for the original model)
One thing that I haven’t specified in my previous post is that I set the top layer of my model as phreatic (the following layers as dependant). Because I do not want my mesh to deform.

Also, I have chosen to treat the “storage change in phreatic top layer where the water table exceeds surface” as confined.
In the original watflow model, “groundwater fow in the elements above the water table is impeded through an ad-hoc exponential relationship between the pressure head p and some artificial "saturation" S:

S(p) = Sr + (1 ¡ Sr)e[sup]?p[/sup]  for p < 0

Where Sr is a residual “saturation” and ? is a parameter controlling the exponential drop-off from full to residual “saturation”.  Partial “saturation” of the elements results in a lower hydraulic conductivity, impeding groundwater flow. This is accomplished by scaling the relative permeability of an element according to its “saturation”: kr = S(p). The conductivity K is now replaced by krK.
Similar to the deforming mesh algorithm, the flow equation is solved with hydraulic head as the state variable. Convergence criteria are imposed for hydraulic head and the "saturation" updates.”

So in the original model I input a residual saturation (0.1 (unitless)) and a fitting parameter (0.2 (unitless)). From my understanding theses parameters would somehow relate to the scaling of sources/sinks by pseudosaturation which is defined through the “residual water depth for unconfined layers” in meters.

As suggested by other posts on the forum, I have tried increasing the value of this parameter.

What I have noticed is that higher the residual water depth value, the less oscillations I have in my errors from one iteration to the other, the error reaches a plateau at a lower value and my head distribution is more realistic. It’s good! (Although the minimum error seems to reach a minimum 0.0003 with a residual water depth of 1 (10 and 100 give similar results to 1). But I remain unclear as to what this “residual water depth for unconfined layers” does exactly, how it affects the rest of my model. Is there a maximum value? How does it relate to the residual saturation?

Thank you again for your help.

Dear Debora,

If you are looking for a mathematical formation, which expresses the non-linearity between saturation and pressure, then I would recommend you to switch to the Richards-equation in FEFLOW. Here you will need to apply an empirical model, e.g. van Genuchten, to describe reduction of conductivity and others behaviors.

About the residual water depth, this could be interpreted as the residual water content. In some Soil Physics books, you will find that saturation can be expressed as the so-called saturation length or saturation depth. If you take an one-meter soil column and you include only the residual water content into it, the column of water inside this container is named residual water depth.

Cheers,

Carlos

Dear Carlos,

It is still ambiguous as to what the “saturation length or depth” means.  I consulted as few soils physics book and I have only come across residual saturation defined as a ratio or a percentage.
In your explanation you refer to a one meter tall soil column, would this mean that feflow calculates the residual saturation as follows:

Sr= residual water depth input by users(m)/1m

In this case the residual water depth would have to be < 1 … however feflow doesn’t limit the input value … inputting 1 or 100 give sensibly the same results.

Thank you for the help,
Debora

[quote author=Blair link=topic=2084.msg4724#msg4724 date=1400116713]
Recharge is assigned using the material property "In/out flow on top/bottom" rather than a BC.
[/quote]
One last question ... how about in 2-D is the recharge also a material property (if so, which one?) or do we apply it as a 2nd type Boundary Condition?

All the help is greatly appreciated!

Debora

2D vertical or 2D horizontal?

I assume you're referring to 2D vertical (cross sectional model) in which case recharge would be assigned through a Fluid Flux BC along the top nodes, rather than the in/outflow material property, which applies to element faces.

My next question would be: is your model solving the standard saturated groundwater flow equations, or Richard's equation for unsaturated/partially saturated conditions? There are difficulties with using fluid flux BCs in 2D unsaturated models as recharge is usually applied at the surface, which is normally unsaturated.

Hi Blair,

My 2-D model is vertical indeed. I am simulating Richard's equation for unsaturated/partially saturated conditions. I did apply the recharge as a Fluid Flux BC. However my model gives heads that are higher than the highest elevation of my model … I am suspecting that it is somehow related to the recharge. I have tried applying a maximum hydraulic head constraint to my Fluid Flux BC but it did not make a difference…

Thanks for helping

Debora

The problem is that when saturation < 1, relative conductivity is reduced according to the Van Genuchten equations (or whatever method you have selected). So the large heads are the result of introducing water mass into a very low (relative) conductivity material.

There are two approaches to manage this:

1) apply recharge BCs just below the modelled water table where s = 1.
2) adjust the Van Genuchten parameters so that Kr is not reduced so severely above the watertable.

I'd suggest you have a read of the help or documentation for unsaturated flow modelling.