- View SourceDear Friend,

I have been compiling selected extracts on groundwater

modelling discussion from various mailing lists for

past few months, as given below. I understand that in

view of random and discontinuous placing, many parts

may not be easily comprehensible. Still I think that

it could be helpful to groundwater modellers to some

extent.

Regards

Kumar

==================================================

C. P. KUMAR

Scientist 'E1'

National Institute of Hydrology

Jal Vigyan Bhawan

Roorkee - 247667 (Uttaranchal)

INDIA

Web Page : http://www.angelfire.com/nh/cpkumar/

==================================================

Kumar Links to Hydrology Resources:

http://www.angelfire.com/nh/cpkumar/hydrology.html

Google Groups in Hydrology:

http://www.angelfire.com/bc/nihhrrc/google.html

Software User Groups in Groundwater Modelling:

http://www.angelfire.com/bc/nihhrrc/groups.html

==================================================

-----------------------------------------------------------------

* In my experience, PCG is faster than LMG for small

models but LMG is faster for large models. Of course,

for small models, speed generally is not so much of an

issue.

* DAMP is used as a multiplier to scale the head

change at each iteration, not RELAX. The DAMP

parameter in the PCG solver package is analogous to

the ACCL parameter in the SIP package. Although I

cannot explain precisely what RELAX controls (since

I'm no mathematician), the typical range for RELAX is

0.97 to 1.0. In my experience, setting RELAX=0.97 can

help to achieve convergence, and may speed up

convergence on a "well-behaved" model, but not always.

* I'm sure you've heard "All models are wrong but some

are useful." If you replaced all the wells in a model

with drains, you would expect to get different

results. How is this any different? You have

replaced one package with a different one that is

based on different assumptions so of course, they get

different results. If you want to know which one is

better, examine the assumptions behind the different

packages and see which assumpptions best match the

conditions at the particular region to be modelled.

If none of the available packages do what you need,

write your own package to do what you need.

* The newer LPF package is the way that most other

models work, i.e., you enter Kz and the inter-cell

conductances are computed internally. The way that

MODFLOW88/96 did it with a pre-computed vertical

conductance always struck me as a bit odd; however,

that has been our happy MODFLOW world for about 20

years. So now we have a new option that is

mathematically more correct, but it is still

disconcerting when the computed heads are different.

I guess you'd have to ask a lawyer about who would

win. My sense is that the differences could never be

explained to a jury. It does present the MODFLOW

modeling community with a problem, though. The way

that you did it with MODFLOW-SURFACT is probably the

way that the USGS should have done it as well - adding

the input of Kz in the BCF Package.

* Input of vertical hydraulic conductivity has its

advantages, however, McDonald and Harbaugh got it

correct in keeping a fixed thickness of a cell for the

vertical direction conductance - for

gravity-segregated vertical equilibrium (GSVE), you

have a transmissivity that is conductivity times

saturated thickness which is valid only for the

horizontal direction and not for the vertical

direction (in actuality, the vertical direction

physics changes from GSVE to unsaturated flow which

would lower the conductivity further and that would be

closer to the "fixed thickness" answers that do not

reduce separation distance between 2 cells vertically

even for consideration of a homogeneous case, but that

aside). Nor does the document state where or under

what conditions the use of a "variable thickness for

vertical conductances" option performs better. That's

why I was asking to see if anyone knew why this was

done in the first place.

* I propose that the node location of a cell that

contains the water table should be located at the

water table instead of the center of the cell as is

standard. Comparisons of this formulation to analytic

calculations numerically implemented in the WATQ code

indicate a much improved solution at shallow depths

compared to the calculations of the current LPF

formulation. This formulation does not, however, take

into account the dynamics of drainage from the

unsaturated region above the water table which can be

a major influence. I want to repeat Richard Winston's

point that "All models are wrong". A model is not

intrinsically good or bad. That judgment depends on

what the model is used for. What may be better in one

situation may not be better in another. You have

highlighted this concept in your discussion of GSVE.

MODFLOW still exists for two contradictory reasons; 1)

as the name implies it is modular and hence easily

extensible to specialized situations. 2) It is widely

known and accepted that the code has accurately solved

groundwater flow problems. It is incorrect to extend

the second point to expect MODFLOW to accurately solve

all groundwater flow problems. Once again required

accuracy like good and bad depends on the problem that

needs to be solved. That MODFLOW produces two results

for the same geometry and boundary conditions means

that it approximates the physics of the flow situation

differently. It proper care is taken, the ability to

change the approximation of the physics of flow is a

good thing because the code can be tailored to

specific problems.

* Yes, having flexibility in averaging schemes is a

good thing. The LPF package provides one more

averaging scheme in the vertical direction. It does a

harmonic average of the conductances, as opposed to

the BCF package's harmonic average of vertical

conductivity times area divided by separation distance

as provided in the MODFLOW documentation for the

various fully three-dimensional, and

quasi-three-dimensional alternatives. However, I am

guessing that people are encountering drastically

different and surprising or unexpected results is

because of the vertical direction treatment.

Specifically, the cell thickness varies as a function

of its saturated thickness in the LPF package for

vertical direction treatment. For a simple example,

consider a homogeneous case with two cells where the

cell below is getting drier. The vertical conductance

quickly rises due to the cell below drying up, to such

levels that drain and dry up the cell above (which

also causes the conductance to increase in a cascading

manner) which would be a very different result than

with use of any of MODFLOW's BCF options. Like I

stated earlier, infact, when things dry up, it should

slow down the conductance due to physics change to

unsaturated flow, and even a gravity flow assumption

will not speed it up.

Yes, there is a dilema in GSVE models for vertical

direction treatment because GSVE by its definition,

applies only to horizontal direction flow - I have

seen this in saltwater-freshwater GSVE models with

multiple aquifers as well. And there too, when you

give it specialized treatment, it is due to certain

physically motivated considerations. Otherwise, the

way MODFLOW BCF did it (or a use of grid-block cell

thickness in the LPF package) is the best you can do

with the VE assumptions being used, since conductance

may be smaller than saturated levels due to

unsaturated zone physics - but you don't know by how

much, unless you simulate unsaturated zone physics.

Using saturated thickness instead of cell thickness in

the vertical direction leads to results that are more

removed from this physics.

* This is one of the most challenging groundwater

modelling projects you can find in real-world

modelling... Ideally, you should use a

saturated-unsaturated code, such as FEFLOW or Modflow

surfact or similar. Such codes are able of simulating

perched water conditions.

If you are using traditional MODFLOW, and assuming

saturated conditions, it is a bit more complicated to

simulate such a scenario. You will probably have some

dry cell issues (and perhaps some convergence

problems), but it is possible to simulate the open pit

crossing multiple layers using drain boundary

conditions.

For the upper layers, where seepage faces may occur,

set the drain elevation slightly above the layer

bottom elevation (say 0.1 m above cell bottom

elevation). If you use Visual MODFLOW as the

interface, that's very easy to do using formulas that

will assign the correct elevation even in irregular

layers. In that case, if the head in the pit wall is

higher than the cell bottom elevation, Modflow's drain

package will remove the outcroppng groundwater. The

conductance value (another input to the drain package)

is tipically a calibration parameter, obtained through

matching drain removal rates with real values observed

at the mine site.

I suggest that you also use a solver that is not too

aggressive, such as PCG2, and that you play a bit with

re-wetting parameters and dumping factors until you

obtain a stable solution. We've done several mine

dewatering modelling projects using this approach and

it definetely works for most cases, as long as you

understand the code and it's limitations.

* Dry cells in a multi-layer model are not necessarily

a bad thing unless they do not accurately reflect the

conditions in the field. Perhaps I am stating the

obvious, but just for the record, dry cells occur in

MODFLOW models when the calculated water table is

below the bottom elevation of the grid cell. Since

MODFLOW only simulates saturated groundwater flow, it

does not represent a head value in the unsaturated

cell, so the cell becomes 'dry' and MODFLOW assigns a

'dummy' head value usually in the range of -1e-30 (or

something else easily recognized as a dry cell).

Although the presence of dry cells in a model is often

a hassle because they create problems with the

convergence of the numerical solution, dry cells do

not necessarily indicate a 'problem' or error with the

model. In your case, the dry cells don't appear to be

causing a problem with the solution convergence (or

atleast, you haven't mentioned this problem) and you

are getting a good mass balance, which is a good sign

but doesn't necessarily indicate 'correctness' of the

solution. If the water levels in the second layer of

your model match closely with the observed field

conditions, then the occurance of dry cells in layer 1

appear to be reasonable.

However, if the dry cells are occurring where you

would otherwise expect saturated conditions (based on

water levels in the field) then you probably need to

revisit your boundary conditions to make sure they are

reasonable (conductance values are often the culprit),

check your initial heads and make sure you are not

starting the iterations with dry cells, and make sure

you have enough recharge flux entering the system to

accommodate flux leaving the system through wells,

rivers, and other boundary conditions. If all this

looks good, then Mark Wilsanac gave some very good

suggestions on how to change the solver and rewetting

package settings in order to deal with dry cells. -

Patrick Delaney

* If the river boundary cells are the only place in

your model where water may leave the system, then one

of the most likely causes of the mass balance problem

is the conductance values you are using for the river

cells. Check your conductance values to make sure

they are within reason. The conductance values can be

reasonably estimated using the following formula:

C = (L x W x Kv)/B

where

L = length of the river in a cell W = width of the

river in a cell Kv = vertical hydraulic conductivity

of the riverbed in a cell B = thickness of the

riverbed in a cell

Most of the popular graphical interfaces for MODFLOW

allow you to enter these 'physical' parameters and

will automatically calculate the conductance values

for each river boundary cell.

Another possible cause for mass balance problems could

be the solver you are using (SIP often produces larger

mass balance than other solvers like the PCG or LMG)

or perhaps the convergence criteria for the solver is

not 'tight enough'.

Just as an afterthought, the fact that you have more

water entering the system from the rivers than leaving

the system through the rivers does not, in itself,

indicate a problem with the water balance. This

indicates you have 'losing rivers' because the water

table adjacent to the river is lower than the river

stage. The excess water entering the system from the

rivers may be leaving through other avenues such as

evapotranspiration, pumping wells, drains, specified

head head boundaries, etc. - Patrick Delaney

* Check out if you dont have any dry well cells, or

any other dry boundary cell. If a cell goes dry during

simulation, it becomes inactive for MODFLOW, and any

boundary assigned to that dry cell will be ignored.

That means that your wells may be inactive (even

though you assigned them to be pumping all the time).

That typically produces a noticeble change in the

budget, but sometimes the heads may not change too

much, unless you look at specific well cells very

closely.

My suggestion is to first look into the budget for

changes in well pumping/injection rates. If the total

rate is different from those you entered into the

model, then you have a dry well cell. Zoom in all well

cells to identify which ones have become dry during

that stress period (look in all model layers).

* With a model set up with cell-to-cell size changes

below 50%, there can be instabilities or failure to

converge. This can be due to other reasons, or cell

size change coupled to highly variable conditions in a

transient solution. If a geologic unit pinches out,

that should be represented primarily by changes to

cell properties, and secondarily by changes in cell

dimension (aspect-ratio rules apply vertically as well

as horizontally). Remember to preserve aspect ratios

of each cell below 10:1, and below 5:1 to be safe. It

is acceptable to have a layer with variable properties

from cell-to-cell. It is rarely acceptable to use

MODFLOW cell dimensions to closely simulate physical

variations in layers when they pinch out.

* I have not used the WHS solver sold by waterloo

hydrogeologic. while i therefore cannot comment on

this solver and its claimed superiority to SIP, i feel

quite sure that SIP has been the best performing

solver for almost all models i have constructed over

the past few years, and never had problems with mass

balance which could not be successfully overcome by

using SIP. i do recall a paper published in Ground

Water years ago when WHS proponents at software

developer Waterloo Hydrogeologic compared their WHS to

SIP and PCG (among other public domain solvers). i

felt that sip had been shortchanged and misrepresented

in their analysis in that they had claimed that sip

yielded a much greater mass balance error than that

obtained when using their whs solver. what they had

failed to do in their analysis, was to accordingly

decrease the HCLOSE with their corresponding decrease

in the acceleration parameter. this is the equivalent

of not closing tightly enough on heads, and therefore,

ending up with an excess mass balance error. this was

an improper use of SIP.

I have found that lowering of the acceleration

parameter, and corresponding decrease in the HCLOSE

will yield a very low (acceptable) to zero mass

balance error, but one must be careful in that the

acceleration parameter is not set so low that the

model converges prematurely with insufficient head

change.

* Check all inputs in the flow model first. Make sure

your flows are realistic. Make sure your K's and

recharge are realistic. Refine your grid and use

dispersivity values that will minimize numerical

dispersion. Make sure your grid cell size complies

with appropriate Peclet criteria (as a rule of thumb,

for finite diferences, Pe should be ~ 2; that means

your grid size should be roughly 2* your dispersivity

value, in the area where, for MOC, this can be

higher). I suggest your read some transport modelling

books that explain this in detail. Be careful with

using constant concentrations, they can act as a

source, but also as sinks (if he concentration in

neihbouring cells is higher than the ones you

specified in the constant concentration node).

* When she is referring to inflow, in plane, and

outflow, I think what she is doing is looking at the

color coding for the velocity vectors or pathlines.

When she sees that there are very few arrows or

pathlines with a green color (indicating 'inplane' or

flat horizontal flow in the layer) she gets concerned

for some reason. Most likely the lack of 'in plane'

flow has very little to do with anything she is

concerned about, it just reflects the fact that she

has slight vertical gradients in the model.

WRT to the vanishing concentrations, you could also

suggest her to make sure she didn't input a default

degradation rate of 0.5 or something very high and

then forgot about it, or the other possibility is she

has assigned a sorption coefficient of 1 (mistaking it

for a retardation coefficient), or perhaps just a very

large sorption coefficient of 0.01 L/mg which would

cause the plume not to move at all, thus giving the

impression that it is disappearing from the model,

when in fact it hasn't moved at all. Another

possibility is that she initially assigned

non-conservative values of sorption and decay

coefficients when she set up the transport model

scenario, and she is now trying to modify these values

in the Setup/Transport Engine dialog instead of in the

Input mode (I've seen this problem a few times) or

Visual MODFLOW.

Finally, the fact that the solution solves in two

iterations seems to indicate that the total run time

of the MT3DMS solution may have been left as the

default of 1 day (perhaps she assumed it automatically

picks up the end time of the flow simulation). She

can change this by selecting Run/MT3DMS/Output Times

and entering the desired Simulation Time.

* The simplest way of defining a starting

concentration is to use the Recharge Coverage and

specify a polygon shape that denotes the contaminant

concentration to be applied to the model.

* What to do when your contaminant plume does not

migrate as expected.

When using MT3D or RT3D for a contaminant transport

simulation, if the plume fails to move as expected (or

does not appear at all in your Visual MODFLOW output),

here are a few common causes and their solutions.

Problem #1: Overlay is not active, or masked by

another overlay

Solution #1: As with Head Equipotentials, Particles,

and even Basemaps and other overlays, you must make

the Concentration overlay active in order to see it.

Click the F9-Overlay button at the bottom of your

screen, and ensure the overlay has been checked off

(to make it active). It is also possible that other

overlays are masking your concentrations. To resolve

this, change the Overlay Order setting to User

Defined, then highlight the Concentration overlay, and

use the arrow buttons to move it up the list.

Problem #2: Simulation Output time

Solution #2: In Visual MODFLOW, the Simulation Output

time can be different for your flow and transport

simulations. Check to see that you have defined the

correct simulation time(s) in the Run Menu, by

selecting MT3D (or RT3D) / Output - Time Steps.

Additionally, in the Output Menu of Visual MODFLOW,

ensure that you have selected the Concentration

overlay as the active overlay, then select the Time

button in the left toolbar, to change output times as

desired.

Problem #3: No concentration input assigned, or not

assigned properly

Solution #3: Ensure that you have correctly defined at

least one cell with an Initial Concentration, or a

transport boundary condition of one of the following

types: Constant Concentration, Recharge Concentration,

Evapotranspiration Concentration, or Point Source.

Please NOTE that Concentration values entered for

Concentration Observation Wells are used for

calibration purposes only; they do not contribute mass

to a transport simulation.

Problem #4: Inadequate flow gradient

Solution #4: Before running your transport simulation,

run just the flow simulation (MODFLOW). Using

Pathlines or Velocity Vectors, ensure that there is a

sufficient flow gradient to cause plume migration.

Problem #5: Incorrect reaction parameters

Solution #5: Check the reaction options in the Main

Menu, by clicking on Setup / Numeric Engines /

Transport. Check that the correct Sorption and

Reaction options have been selected.

First order decay rates (dissolved and sorbed phases)

may have a tremendous impact on the plume mass, since

during the simulation, some contaminant mass may be

removed from the model. If the decay rate is too high,

your plume will show very small concentrations, and/or

will not be visible at all at later simulation times.

Decay rates (due to biodegradation or other

mass-removal processes) can be taken from literature

values, but this should be done with caution, as this

parameter is highly site-specific for most compounds.

A good reference on decay rates used in natural

attenuation studies can be found in:

Wiedemeier et al., 1999: Technical Protocol for

Implementing the Intrinsic Remediation with Long-Term

Monitoring Option for Natural Attenuation of Dissolved

Phase Fuel Contamination in Ground Water. Air Force

Center for Environmental Excellence, Brooks, AFB.

This document can be downloaded from:

http://www.afcee.brooks.af.mil/products/techtrans/monitorednaturalattenuation/protocols.asp

The decay rate (l, or sometimes called Kd, with units

of 1/day), is typically obtained from half-life

values, converted into appropriate units using the

following relationship:

l = ln(2)/t1/2

Where t1/2 = half-life of the compound

In addition, check that your Kd value is correctly

defined in the Species Parameters tab. A very large Kd

value will result in an extremely high retardation

factor, which can result in lack of contaminant

movement through your model.

For typical organic compounds (such as TCE, DCE, PCE,

BTEX, etc.) where linear, reversible sorption can be

assumed, retardation will be calculated using the

following formula:

Retardation = 1+ (Bulk Density/porosity)*(Kd)

Where Kd = Partition coefficient

For hydrophobic organics, Kd can be determined using

the relationship below:

Kd = Koc*foc

Koc - octanol-carbon coefficient

foc - organic carbon fraction in the aquifer

A more detailed overview on Kd's can be found in:

US-EPA, 1999. Understanding variation in Partition

Coefficient, Kd, values. EPA

402-R-99-004A&B. US-EPA Office of Air and Radiation.

This document can be downloaded from:

http://www.epa.gov/radiation/cleanup/partition.htm

* From the Modflow standpoint, water that exits the GW

system through of a drain cell is accounted for in the

budget, but is otherwise not simulated. The situation

for river cells is similar--Modflow does not simulate

any flow in the river, it only simulates exhange

between the GW system and each individual river cell.

So there is no way to route water from a drain cell to

a river. The Streamflow Routing (STR) Package, on the

other hand, does simulate flow in surface-water

streams, accounting for flow rates in the stream and

calculating stream stage, which is then used as the

external head in the GW/SW exchange calculation. If

you want to simulate exchange between drains and

rivers, you could use the STR Package for both types

of features.

* How to import modflow model files to gwvistas 4?

Select File/New and click the MODFLOW button. On the

next dialog, click browse to go find the name file

(*.nam). If this is an older MODFLOW88 dataset, then

browse to find the *.bas file. Some datasets

(especially those that were created without the use of

a preprocessor) can be troublesome.

* I want to input separately 2 sets of recharge into a

modflow model: (1) Effective recharge from rainfall,

and (2) urban leakage (same format as the (1)), which

will overlap. Since both components are fairly complex

and over a long time (1200 Stress Periods), I would

like to keep them as independent as possible (avoiding

the data compilation option, if you see what I mean).

It is possible to do what you want using just the

Recharge package. Here is how you do it.

Define each recharge distribution for each stress

period using parameters. You can use one set of

parameters for the rainfall and another for the urban

leakage. With each parameter, you can use multiplier

arrays and zone arrays to define the spatial

distribution of the recharge. Then, for each stress

period, use the parameters for that stress period to

apply recharge from both rainfall and urban leakage.

MODFLOW-2000 will add the contributions from all the

parameters you use in a stress period to determine the

total recharge.

You should be able to do this with any reasonably

up-to-date GUI for MODFLOW-2000. If you are using a

GUI and you can't figure out how to do this, contact

the GUI developer for assistance. If you aren't using

a GUI, you can check the input format in the Online

Guide for MODFLOW-2000 at

http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/guide.html

or in the original MODFLOW documentation as modified

in "Time-varying-parameters.pdf" (distributed with

MODFLOW-2000).

* Since I know you are using Modflow VKD (based on

MF96, not MF2K) and that the wells package is almost

certainly being used for many abstraction wells, I

suggest 2 routes:

1. A FORTRAN utility to combine two recharge files - a

days work? 2. To use a MF package that you aren't

already using - but with some manipulation e.g. the

rivers package. By setting the stage of the river as

above any elevation in the model (1000m?) and the

conductance equal to the quantity you want to recharge

for that SP - with a 1 m difference between stage and

bottom elevation you will get the appropriate amount

leaking (recharge) to each cell.

You can achieve the same result with the drains

package by specifying 2 drains in each cell. The trick

is to have one with a positive conductance and the

other negative, but the same value. Make the

difference in elevation 1m and the magnitude of the

two conductance values the quantity you want to

leak/abstract. If the elevation is outside the range

of variation in the model, you end up with a constant

recharge/abstraction from the model.

For composite/trick boundary conditions of this type I

suggest drawing figures for them like those in the

MF88 book, it will help to understand what I'm trying

to say.

* Design the shoreline?

Depending if you are using FEMWATER or MODFLOW, it is

handled differently. With FEMWATER, you want the

boundary of your 3D mesh to follow the shoreline. With

MODFLOW, you are working with a 3D grid. So, you need

to mark the cells that contain the shoreline with a

fixed head boundary condition, and the cells outside

this area will be disabled.

If you are attempting to model the interaction with

the sea bed, then you will need to set the top

elevations of the mesh or grid to match that of the

sea bed. However, generally you end your model at the

shoreline.

* MT3D Performance

What you need most depends on where your bottle-neck

occurs. You need to monitor CPU and RAM usage on the

computer during the MT3D run. If you run Windows 2000

or XP, you can right click (for a right-handed

pointing device) on the Taskbar (grey bar at bottom of

screen), then select Task Manager in the menu pop-up.

Once Windows Task Manager opens, select the tab

labeled "Performance". "Performance" shows CPU usage

graphically and memory usage in numeric form. If the

"CPU Usage History" (top graph) stays maxed at 100%,

then a CPU speed increase should help. In parallel,

observe the "Physical Memory" values. If the

"Available" value becomes a small fraction of the

"Total", AND the "Commit Charge, Total" exceeds the

"Physical Memory, Total" then you are likely to have a

RAM-limited condition. In this case, more RAM would

help. With 300,000 cells, consider 1 Gbyte or more. If

both things are happening (unlikely), then upgrading

CPU and RAM should help. Be aware that upgrading one

component sometimes results in the other becoming the

limiting condition. That is why it is "unlikely" both

limits are happening simultaneously. Usually one tops

out before the other. If you run Win98... consider

upgrading the OS.

* There are various things that you have to do to get

Surfact running in Vista 3:

1. You can only have a limited number of packages.

For example, if you are doing contaminant transport

modelling then you have to turn some other pacakges

off, otherwise Surfact will not run e.g. you can turn

the cell by cell flow packages off (just type 0 into

the packages number in Model - Modflow - Packages).

Another pace you can turn packages off is in the

output control: Turn off the drawdown unit.

2. You have to modify some things in the Model -

Modflow options and some in the Modflow - Surfact

options. These are:

Model - Modflow - Packages: change modflow version to

surfact; change solver to PCG4; and untick

automatically reset package units

Model - Modflow - BCF package: tick the box to use

the BCF4 package

Model - Surfact - Packages: add a unit number (one

that you have not yet used) and tick the box for BCF4;

add the unit number to PCG4 and tick the box (use same

unit number as shown in Model - Modflow- Packages);

add a unit number and tick box to RSF4 is you want the

model to simulate surface seepages (and remeber to

turn this unit on in Model - Surfact - RSF4 package).

Model - Paths to models: set modflowwin32 option to

DO NOT USE; put in correct path for surfact in the

modflow box.

I can't think of anything else right now - I always

forget one thing and spend a couple of hours wondering

why surfact won't run. If you have any difficulties

it would be helpful to know if the surfact Dos window

displays and if so, what it says in the window.

* The most common problem when using SURFACT is

setting up the program file. Select Model/Paths to

Models. Change the MODFLOWwin32 Option at the top to

"do not use". That tells Vistas that you are not

using one of our model DLLs. Then change the MODFLOW

program from MFWIN32.DLL to whatever your SURFACT

program is (usually it is msft.exe).

* I believe that if you are using the latest version

of SURFACT (V2.2) that problem 1. below is no longer

an issue. It used to be that SURFACT could only open

a limited number of files but that is no longer the

case. Also, if you upgrade to Vistas Version 4 and

SURFACT V2.2 the link between the two programs is

easier to establish. The same holds true for all of

the other versions of MODFLOW that Vistas supports

(MODFLOW88, MODFLOW96 in double precision,

MODFLOW2000, MODFLOWT, MODFLOW-SURFACT, SEAWAT, and

SEAWAT2000).

* VS2DI or VS2DT, being a Richards equation-based

models have certain PROs and CONs when used for

groundwater recharge estimation:

Pros: the Richards equation models are the most

theoretically based and allow representation of the

flow processes in porous mediums more realistically

than the water-balance models. They have been well

tested against field and laboratory experiments and

proved to provide good correlation with observed

results.

Cons: the major complicity of the Richards equation,

comparing to the saturated flow models, is that its

coefficients are non-linear functions of the pressure

potential. To approximate these functions, three

different algebraic equations are usually used (named

after their authors): Brooks and Corey's, Gardner's,

and van Genuchten's. Richards equation can be solved

with a very limited number of boundary conditions. The

condition for water at the boundary can be specified

either as the flux of water across the boundary, the

head at the boundary or as a combination of specified

head and specified flux. The Richards equation with

these boundary conditions is a nonlinear partial

difference equation that has no close-form or

analytical solution. To solve the Richards equation,

a finite-difference method of numerical approximation

is applied in VS2DI. Searching for a best numerical

method for the Richards equation became a special

field in hydrologic science, particularly in 1980's,

however, almost all implemented approaches cannot

assure that the model will unconditionally converge

and simulation results will be obtained. It is an art

to get the model to work correctly (produce results

for the whole simulation period with a good water

balance) when you have varying flow boundary

conditions in your profile. Our experience also proved

that applications of Richards equation-based models to

highly heterogeneous soils with variable hydraulic

properties can be extremely difficult to execute and

time consuming.

These days there more and more examples when less

sophisticated in flow equation but much more advanced

in terms of weather/plant effects water-balance model

HELP is used to estimate groundwater recharge.

* I will suggest you keep the north and south boundary

as variable head boundary. This I am suggesting you

treating there is no river, lake etc. If these

features are present on north or south side of model

domain, then you have look for another type of

boundary. Regarding limited observation well data and

that too is limited to layer1, then it will not be

possible to do any realistic modelling. If there is no

scope of generating additional data, i would like to

suggest that you confine your modelling activity to

first layer only.

* What makes a good correlation depends on what your

expectations are. If you are conducting exploratory

research, maybe 0.2 isn�t bad. If you are testing a

known process having some natural variability, 0.6

might be acceptable. But if you are calibrating

instruments, maybe 0.9 isn�t good enough.

If you don�t have any clear expectations, how do you

tell if a correlation is good? The answer comes in

three parts.

1. Value - Square the correlation coefficient (called

the coefficient of determination or just R-square).

R-square is the proportion of variation in the

dependent variable (y) that can be accounted for by

the independent variable (x). You might be able to

decide how good your correlation is from a gut feel

for how much of the variability you wanted your model

to account for.

2. Significance - Every calculated correlation

coefficient is an estimate. The �real� value may be

somewhat more or somewhat less. You can conduct a

statistical test to determine if the correlation you

calculated is different from zero (or some other

number if it�s relevant). The larger the calculated

correlation and the greater the number of samples, the

more likely the correlation will be significantly

different from zero. For example, a correlation of

0.59 (R-square of 0.35) would be significantly greater

than zero based on about 25 samples, but a correlation

of 0.01 wouldn�t be significant with 250 samples.

3. Plots - You should always plot the data used to

calculate a correlation to ensure that the coefficient

adequately represents the relationship. Specifically,

the data should be linearly related and free of

outliers. There are a few other things to look for

(e.g., hidden trends, autocorrelation) that I won�t go

into here.

So, the reviewer who said that �the "R square" value

of 0.35 has no significance and is just as good or as

bad as say 0.01� would only be correct if she looked

at a statistical test of whether the value was

significantly different from zero. R square does not

need to be more than any particular value to draw a

comparison. Remember, too, that �no relationship� may

also be an important finding.

* If you are considering GUI's for your project, you

may want to look at Visual MODFLOW Pro. In

particular, if cell splitting, or grid refinement is

an issue for you. Visual MODFLOW includes two

distinct features that help users with this.

1) Cell refinement tool

2) Grid smoothing tool

1) Cell refinement/splitting tool: You'll need no

more than a few clicks to edit your grid in a very

flexible manner. To split your cells, just click on

'Refine by ...", and specify how many times you want

the cell to be split into (2, 3, 4...). Then, select

the area you need to refine with your mouse.You can

apply this to one row or column, or all cells in the

domain - whatever you wish. The same applies for

splitting layers and coarsening the grid.

2) Grid smoothing tool: After grid refinement, you may

want to use the grid smoothing routine to produce a

smooth transition between coarser areas to highly

discretized ones. Again, just a few clicks. This is

actually a perfect tool to help you better understand

the effects grid size and refinement will have on your

model results.

Another feature that may be useful for you is Visual

MODFLOW Pro's batch capabilities. It will allow you

to prepare all data input files in different model

projects and run them in 'batch mode'.

* The interpretation method is different according to

the geology. You have methods for porous media, others

for fissured aquifer, others for double porosity

media. Methods also vary according to the type of

aquifer (confined or unconfined; semi-confined with a

leaky aquifer...) and corrections can be made

according to the type of well (partially penetrtaing

well; pumping well with head loss; skin effect and so

on). Finally, boundary effects (barrier or recharge)

may affect significantly your interpretation.

In order to remain pragmatic I would suggest you to

start with the simplest method : the Theis formula,

(or even the Jacob approximation if your pumping time

is big enough) and check wether you get a good fit or

not with realistic parameters. Even a fissured

unconfined aquifer can be interpreted with Theis when

your radius of influence is wide so that the fissure

network can be assimilated to an homogeneous

equivallent porous media and the depletion of the

water level is moderate compared to the thickness of

the aquifer. This should give you a rough estimate of

your transmissivity, which is in many cases,

sufficient.

If you need more accuracy you can use specific

softwares. You will find on the net many softwares

proposing many different methods. You don't need Pest

as they have their own fitting engine. AQTESOLV for

instance proposes methods with well storage effect.

However, as soon as a software proposes its own

adjustement, it's becoming very difficult to know

what's you're doing and you end up playing a video

game. So be carefull, especially if you lack

experience in this domain, not to use complex methods

with many parameters which will end in a perfect fit

but with irrealistic results. In my company (French

geological survey)we use ISAPE, a software which is no

more commercialised, but whose philosophy is to let

the user propose realistic values and where you can

add or remove well effects or boundary influence. I

wish this approach were more widespread.

* I think your scenario will represent two different

geologic materials with peculiar intrinsic properties.

Hence, I suppose this will require two separate

simulations. To the best of my knowledge, MODFLOW 2000

only allow one K matrix values per layer per

simulation and this is incorporated within the flow

package. However, if the introduction of the permeable

barrier is controlled along defined specific paths or

unique relatively smaller area(s), then you can use

one simulation with two stress periods and incorporate

Horizontal Flow Barrier Package in the second period

to account for the permeable barrier distribution.

This is not suppose to be of much problem especially

if you incorporate your modflow with GIS package like

GRASS GIS.

* Sophisticated modeling programs like Visual Modflow

Pro, or any of a number of similar packages, will run

properly with a variety of inputs that may not be

physically realistic - and thus provide garbage for

output. With that in mind, here are some thoughts on

your questions:

1) You may divide a single aquifer into more than one

layer if there is some reason to look at vertical

stratification within the aquifer. For example, a

thick aquifer may have different hydraulic or

transport properties in different vertical zones. You

may also divide a single aquifer into several layers

if you want to examine the vertical flow patterns due

to pumping from different zones within the aquifer.

2) The top of Layer 1 could be the ground surface,

the water table (not as common because it can move

vertically), or the top of the uppermost aquifer or

aquitard of interest depending on the physical system

being modeled.

3) You'll need to assign boundary conditions no

matter what the size of the physical system is. If

you don't have natural physical boundaries, one

strategy is to assign boundary conditions at a

distance far enough away from the area of interest to

minimize the effects of boundary assumptions on the

solution.

* It is appropriate to use a single layer if you feel

that it is reasonable to ignore vertical flow in your

model. Assuming that you would be using the Layer

Property Flow (LPF) Package of Modflow-2000, and if

the aquifer represented by layer 1 is not overlain by

a confining unit, then specifying the top of layer 1

as ground surface is reasonable. The reason that you

would be better off this way than specifying the water

table elevation as the layer top is that the LPF

Package only provides two options with respect to the

confined/unconfined status of a layer: Confined or

Convertible. If the head in a cell goes above the

specified top of a convertible layer, the cell is

treated as confined. Every model needs to have at

least one fixed head, either as a constant-head cell

or as a fixed external head (such as the stage at a

river cell), or the solver will not be able to reach a

solution. You may need to expand the domain of your

problem to encompass a fixed head.

* 1- Model design depends greatly on your project

goals. For example, you could have one aquifer, but an

overlying layer (and/or underlying layer) as well. Or,

you could be modelling just 1 aquifer, but split it up

into several layers in VMOD to provide more detailed

flow information (i.e. if the aquifer is 100 feet

thick, having just 1 layer will not provide very

detailed output. You can subdivide the 1 "real-world"

layer into 10 "model" layers, all with the same

hydrological properties, in order to obtain more

detailed output). Or, you can simply design a 1-layer

model if that is all your project requires.

2- In Visual MODFLOW, Layer 1 represents the uppermost

layer of your profile, and the top of Layer 1 is

therefore the ground surface. Layers in Visual MODFLOW

represent soil layers (or conceptual soil

layers)....the location of groundwater in these layers

is determined by your model inputs, and the resulting

calculations of the flow model.

3- The assignment of boundary conditions is a question

that often goes beyond the scope of Technical Support,

and enters the realm of Extended Modeling Support.

However, to provide you with some guidelines:

-If you have no "obvious" water inputs (streams,

rivers, lakes, injection wells, recharge from

rainfall, etc.), you can try assigning an upgradient

equipotential line, and downgradient

equipotentialline, using Constant Head Cells at the

upgradient and downgradient edges of your model, to

represent the regional water table (other boundaries

may be more appropriate, depending on your specific

model domain and objectives). You can then add

additional inputs (extraction wells, contaminant

sources, etc.) to the model region, which will be

influenced by your regional flow gradient.

4- I would recommend taking some time to run through

the Tutorials included with your Visual MODFLOW

software (the PDF instruction files are located on

your installation CD-Rom, and the files are

automatically copied to the Tutorial subfolder in your

Visual MODFLOW program folder). These tutorials are

designed to teach you how to work with your software,

and how to design a model, import data for model

inputs, run the various packages, calibrate your

model, examine your output in 2D and 3D, and export

data.

* You can use the field interpolator program in the

PMWIN package. Using the field interpolator you can

interpoltae the field data in txt file to the grid

point of your model. On the other hand, you can import

any dxf file which illustrate the changes in hydraulic

properties to help you in the graphical interface.

* Upwind discretization always has an artificial

diffusion term no mater how refined your grid is (even

thought it diminishes with grid refinement). Central

discretization always produces an artificial

dispersion term that causes "wigles" however you can't

control it with grid refinement. For almost any Pellet

number the upwind scheme will produce positive

coefficients for your coefficient matrix, while that

isn't true for central differences. However central

differences has "second order accuracy" while upwind

only has first order. You can visualize numerical

diffusion and dispersion mathematically by deriving

the "equivalent or modified equation" see

http://widget.ecn.purdue.edu/~jmurthy/me608/main.pdf

for some cool notes on this (as good as any book I've

read).

Physically you can visualize numerical diffusion by

imagining a 1D discretization of a domain | : | : | :

| where | are faces and : is the center of the cells.

If your initial property is 1 on a given cell say i

and 0 on all others it would look like this:

| 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0

So using upwind with a courant = v dt / dx of say 0.5

and a velocity from left to right after the first time

step the accurate solution would be something like

(zooming on cell i):

| 0 { 0.5 | 0.5 } 0 |

Note that I have created a new cell delimited by the

{} brackets this "new cell" that doesn't exist in our

domain is exactly half of cell i and alf of cell i+1.

The property's value will be 0.5 inside this new cell

and 0 outside it (meaning that cell i and i+1 would

have half with 0.5 and half with 0). This would be the

accurate solution because a courant of 0.5 makes the

property shift half a cell every time step. However we

are bounded to our initial discretization so the

solution we will obtain will be:

| 0 | 0.5 | 0.5 |

so you can see than instead of having a 0.5 of

property in a volume the size of a cell we have 0.5 in

a volume the size of 2 cell thus lowering the

concentration. This happened because the definition of

concentration is C = lim vol ->0 DM/DVOL and we

implemented it as DM/DVOL. If we would diminish the

size of our grid so that currant = 1 we could obtain

the exact solution. So the smaller your grid is the

more accurate your solution will be. This is the

physically correct way of solving the diffusion

problem another way is to use higher order schemes

witch usually maintain the sharpness of the gradients

but usually remove mass from all the wrong places

(even thought they conserve it). Looking at the same

example a central differences discretization would

derive something like:

| -0.25 | 0.75 | 0.25 |

Cell i-1 will produce a negative property because the

concentration al the face between i-1 and i is 0.5

cell i +1 will obtain the accurate concentration of

0.25 (0.5 in half the volume of a cell). But we are

producing the famous "wigles" close to the properties

gradient. There are smarter second order schemes like

qwick or TVD schemes.

* In visual MODFLOW there is provision for importing

MODFLOW files (*.bas files). In PMWin there is

provision for conversion of MODFLOW88/96 (*.nam). I

have tried to import *.bas file in Visual MODFLOW and

found that model data is not correctly coming to

model. These options may be available in GMS but I am

unable to locate the exact sub-menu. Number of options

is available for opening different format files In the

FILE menu. But It does not permit opening of *.bas or

*.nam files. This option must be in different menu.

-----------------------------------------------------------------

__________________________________

Do you Yahoo!?

Dress up your holiday email, Hollywood style. Learn more.

http://celebrity.mail.yahoo.com