Jonathan,
'rech' represents a simple source applied only to the top-most face (or cells abutting the top most face) of the model. Its rainfall that recharges water in the subsurface. Probably the best way to account for this would be with a FixedGrad boundary condition along these faces, but I'm under the impression that this boundary condition doesn't work in FiPy, at least not for all conditions. Typically, this data is available with units of length per time (eg. inches per year). My understanding is the source terms in FiPy need to be expressed in units of per time, at least for what I'm doing, and to convert it to an appropriate value I need to multiply this rate by the area of the face it crosses and then normalize by the cell volume. This gets a little messy as one goes from a nicely structured and regular grid to one with variable cell and face sizes This is the reason I resorted to my convoluted approach and didn't just put source term in all of the upper-most cells. It would be great if there is a simpler solution for a fixed flux boundary. Just before I got your message I was looking at this again and saw I also did some other confusing things in the script. I'll rework the script and see if I can't make it a little more readable and see if simplifying this even more allows me to calculate mass balances within the model. Thanks for taking the time to look at this. Andy
