(long) question on phreeqc

[Chris.vanUffelen@xxxxxxxxxxxxxxxxxx]

Dear Mr Parkhurst,

I am a Dutch soil scientist working on my PhD-thesis. 
My special interest is the formation of sodic soils in
Hungary.  I use phreeqc a lot to analyze my data, and
I really appreciate its capabilities.  I have followed
the phreeqc course of Tony Appelo in Amsterdam.

Presently I want to use phreeqc to do calculations on
my concepts of soil formation.  I'm afraid I got stuck
with phreeqc's transport capabilities.

To clearify my problem, I think it is best to describe
how I experience the present TRANSPORT module in
phreeqc.  My concept of transport in phreeqc is as it
behaves like a freight train in which wagons carry the
solutions along the flow line.  Along the flow line,
there are a number of platforms which serve as the
precipitates and exchangeers.  The idea is depicted by
the following (cells 1-5 belong to the column, cell 0
contains SOLUTION 0, which contains the initial
solution, furthermore, `o' implies the original
solution, `n' the new solution, and `X' the front
between `o' and `n'):

cell_no:    0 1 2 3 4 5
step 0:     nXo o o o o
step 1:     n nXo o o o
step 2:     n n nXo o o

As the transport process proceeds, wagons with either
`n' or `o' move along the different plantforms,
equilibrate, and diffusion is taken into account by
mixing the contents of adjoining wagons (here the idea
of a freight train fails a little).  This model
produces good results in case of saturated flow: the
different solutions `n' and `o' are carried in their
wagons (and the front `X' may be regarded as the
engine of the train).  The wagons before and behind
the engine will proceed and there is nothing against
nearly empty wagons (in the phreeqc examples I have
seen, waggons tend to be filled, or, flow is always
saturated flow).  I can depict this in the following
way.  Instead of the symbols `n' and `o' I'll write
down the water saturation in the cell (which for this
simple case is trivial of course, for the water
saturation equals 1 in case of saturated flow).

cell_no:   1   2   3   4   5
step 0:    1   1   1   1   1
step 1:    1 X 1   1   1   1
step 2:    1   1 X 1   1   1
etc.

However, this concept produces unwanted effects in a
sub-saturated soil which becomes saturated as a water
front infiltrates the soil.  For the type of work I do,
there is presently nothing against a mixed cell concept
(I mean, I do not presently need Darcy equation with
K(h) relations that relate the conductivity to the water
suction in the soil).  However, the TRANSPORT model
needs a (little?) modification.  What I need is a system
in which water before the front (so initially present
water) accumulates to saturation 1 before it starts
moving.  For example, a column in which the initial
moisture saturation is a half, the column develops in
the following way when in each step a volume of water
equal to half the saturation in cell is added:

cell_no:      1   2   3   4   5
step 0:      .5  .5  .5  .5  .5
step 1:      1 X .5  .5  .5  .5
step 2:      1   1 X 1   .5  .5
step 3:      1   1   1 X 1   1

This has the consequence that the solution in
cell 4 starts moving in step 3, and that only
after step 3 water will leave cell 5.  It is
in principle also possible to deal with a
slighly more complicated example in which the
initial column has different saturation
values:

cell_no:     1  2  3  4  5
step 0:     .2 .3 .1 .5 .2
step 1(+.8)  1X.3 .1 .5 .2
step 2(+.6)  1  1X 1 .5 .2
step 3(+.3)  1  1  1X 1  1

One advantage of such an approach is that
there is no chance that a highly
concentrated solution at a certain depth
starts moving too early and producing
unwanted effects, or disappears from the
system even though it is not reached by
the infiltration front.  There is,
however, an even more important advantage
to such an approach.  It allows me to
simulate capillary rise due to which
considerable amounts of salt may
accumulate at the soil surface, because
all the water which enters the soil at
the bottom boundary eventually
evapotranspirates somewhere in the soil
column.  I think, I can achieve such a
behaviour by means of a KINETICS and
RATES couple which extracts part of the
water from the upper cells, near to the
soil surface.

I hope I succeeded in describing the
requirement I think I need to fully
exploit phreeqc chemical capabilities. 
Now, here are the questions I would like
to ask:

1  Is such an approach reasonable? I'm afraid
   I'm not a (numerical) mathematician, so I
   find it difficult to judge.  There is a
   problem I foresee, which has to do with the
   additional mixing of solutions in adjoining
   waggons to allow for the water accumulation
   at the front between sub-saturation and
   saturation.  This might interfere with the
   control on hydrodynamic dispersion.  Maybe
   this can be solved by the use of finer grids.
2  If it is a reasonable approach, is it possible
   to implement it in phreeqc?
3  If that is the case, can it be done by
   modification of the present TRANSPORT
   module, or is it better to define a new
   module, e.g.~SOIL_COLUMN, and to what
   extend do you think it can make use of the
   code already present in TRANSPORT.
4  Are you aware of any such modifications, or
   modifications that take into account the aforementioned
   K(h) relationships?

I am also not a programmer, although I have experience
with C.  If it helps and if the possibilities are there
I am happy to do the programming, as I think phreeqc is
a good programme and the way it is offered to the
(scientific) community serves as an example for the way
scientific work should be treated.  If I can contribute
to that, I would gladly do so.

Thanks for your time,

Chris van Uffelen

Sub-dept of Soil Quality
Dept of Environmental Sciences
Wageningen University
The Netherlands