PHAST for multi layer aquifer system

[rapheul]

Hi, there,

I am starter of PHAST, hope my questions are not too basic.
1. I have a 4 aquifers system as well as aquitards in between. Both unconfined aquifer and confined aquifers. Should I use different prisms media to build it up ? Accordingly, with Free Surface BC true, is phast refering to the top layer automatically or? I did not see any detailed definition of the free surface domain.
2. I would also do the transport modeling. And of course each aquifer has different chemical properties. For prism media definition, there seems to be no solution available for definition. How should I set initial chemical conditions for aquifers? But I notice for any BCs there could be a solution assigned. What would you suggest if I do such modeling while I have hydraulic links between aquifers? besides, different aquifer has various BCs, how should I set up different BCs for the system?

Thanks in advance~


 

[dlparkhurst]

PHAST has only two options: a free surface, or not a free surface; it applies to the entire domain.

The PHAST domain is a full 3D volume of space. Unlike Modflow, it is not a series of connected layers. You will need to define zones (possibly prisms) that provide the media properties of the aquifers and the confining layers. PHAST is a little funky in the way properties are distributed. MEDIA properties are distributed by "element", which is the space defined by 8 nodes. Properties are averaged at each node. To be sure you have the low conductivity of a confining layer, you should define a grid that has two vertical elements with the confining unit property.

Initial chemical conditions are defined in CHEMISTRY_IC. It is there you would define solution compositions that would apply to each aquifer and confining unit. It accepts zone input where any node within a zone will be assigned the solution you specify. It is here that you also define all of the other reactions (equilibrium phases, surfaces, exchange, kinetic) that apply in the various parts of the domain.

Flux and leaky boundary conditions are also defined by zone, except that only exterior faces of cells within a zone get the boundary condition. So you can probably use the same prisms that define each aquifer to define the boundary conditions for each aquifer.

Hopefully you are using Phast4Windows to set up your model. It is a qood way to see the various input options and to visualize the zones that you have defined. The current version also allows you to run the model with parallel processing.

[rapheul]

Hi David,

Thanks for your quick answer!
Yes, I am using Phast4windows now. It is much easier to visualize the model. And true I see that phast does not have a connected multi layers. I am sorry but I still don't quite understand the first question. What is "a gird has two vertical elements for confining layer"?  Do you refer to define another single prism vertically between other prisms(aquifers) as confining layer and then assign its K with a very low value?  Could you please explain it more? Thanks a lot.

For boundary conditions, i wonder how I should define if the aquifer(confined) has a specific water head BC? because for my case, i have no data for flux but some wells with obs. heads. I saw in the tutorial there were clear BCs shown in the pics. I'd also like to show them out after definition. How to realize that(Fig 6.5 in tutorial)? It will be also helpful to check whether you defined a correct BC for example the leaky areas.

With best regards and thanks again~

[dlparkhurst]

Please look in the manual for the description of spatial properties. Media properties are defined as the volume weighted average of the elements that meet at a node.

Media properties can be visualized from the HDF file. Run zero or one time step and use Model Viewer to look at the media properties.

You can define specified head boundary conditions.

The boundary condition zones are shown in P4W and you can visualize the nodes that have each type of boundary condition in Model Viewer.

[rapheul]

OK, I will look for those details. many thanks~
Regards

[rapheul]

Hi David,
Sorry to mention this again here. Do you have an example for this kind of multi-layer aquifer model? In the manual, there is no such exercises.  I did it as you told me last time to define different prisms but when I looked at the media property in Model Viewer, there was always a transition zone between, the interpolation of the property values(K) was done automatically but actually I assigned only two constant values for the two prisms. If it should not be connected layers, how can I define those mult-layers system with specific elevations that are connected in real world? It would be better if there is an example p4w that has similar system.
Really sorry for this, but I still did not work out a good way to solve my problem. Thanks in advance.

[dlparkhurst]

As I said earlier (and as explained in the manual), media properties, such as X hydraulic conductivity, are averaged at each node from the definitions for the surrounding elements. Thus, if there is a transition in a property, there will not be a direct change from one value to another. There will always be a node which contains an averaged value.

This is not a problem if you have fine resolution, there will be one averaged value at the transition in a given coordinate direction. This is modeling, an approximation of an aquifer. In any case, you have very limited understanding of the actual hydraulic conductivity field for any aquifer.

Also, if the layers are horizontal or vertical, for a given grid, it is not difficult to ensure that your aquitard is thick enough to have at least one node that has a specified low hydraulic conductivity.

However, if you have a coarse resolution and sloping boundaries, there can easily be problems. Suppose you have two aquifers separated by an aquitard. It can easily be that you will define a zone for the aquitard that contains none or only one element in a coordinate direction, in which case, you may end up with a node that you expected to have the value of the aquitard, but has either the aquifer value or an average of the aquifer and aquitard. So you must be careful when you define these zones.

There is an option to ensure that at least one node has the value defined for the aquitard. It is done with the -shell option of a zone definition for MEDIA. -shell is used to define the properties on the perimeter (surfaces) of a zone (not the interior). In the example attached, a zone with a sloping plane for the top is used. The X hydraulic conductivity is defined to be 0.01 in the "shell" compared to 1.0 for the aquifer. The shell is defined to be 20 m thick in the Z direction. If the grid has 101 nodes in X and Z, then the aquitard is defined by about 20 nodes with a 1-node transition above and below. If the grid has few nodes, say 11 or less in Z, then the zone is defined to contain at least 2 elements in the vertical direction.


(The flow system has a specified head in the upper nodes X=0, Z=100, H=100; and specified head at X=100, Z=80,H=80. It is a steady-state, free-surface calculation. )

[rapheul]

Thanks for your kind explanation again, sorry for my bad understanding of the whole question. until now have I understood the definition of media property. I read it in the manual but not sure whether the distribution should be like that or not.
And by the way, another question.  what would you suggest if I try to set up mining spaces which I will use later for the water flooding question as well as geochemical evolution? means that the mining layer was empty first and then after closure it was flooded with water.  both the aquifers below and above will be the flooding sources. is there a proper definition for this?
Thanks so much for your kind help.
Regards

[dlparkhurst]

The same logic applies to more permeable zones as to my descriptions for less permeable zones. You can define a mined level with high hydraulic conductivity, but you still need to worry about the averaging of properties that PHAST does.

You did not say what kind of mines. But let's assume it is a mine where pyrite can oxidize and cause acid mine water. It may be difficult to account for the mined zones for a couple of reasons: (1) the stopes may have accumulated oxidation products over the mining years that will dissolve with the first flooding of water, and (2) oxidation may occur as long as the stopes are open to the atmosphere, but will cease when the stope is completely flooded.

For (1) the first flush of water may be quite concentrated and acidic. You could use a flux boundary to add water (hopefully you have some samples to work from) for the period of time a stope is at the water table. Alternatively, maybe add limited quantities of oxidation products in the equilibrium phases that will dissolve when the mines fill.

For (2), maybe you can ignore the additional accumulation of oxidized products while the mine is filling, assuming the accumulated oxidation products are the dominating influence. Once the mines are full, very little additional oxidation will occur. Changing the set of reactants in the middle of the run (removing the fixed partial pressure of O2(g) for instance) would be difficult. Again, maybe a limited amount of O2 in the equilibrium phases would simulate the presence for a period of time that ends when all of the O2 is consumed.

[rapheul]

Thanks for your advices. Actually I have the same idea with you, also plan to define a high permeable zone as stope. Then is that correct if I define the initial condition of the zone be totally "dry" (with a 1st kind BC lower head than the bottom elevation? )and then later on with flux BC or well injections to add water? Or do you have better suggestions to make it from"dry" to "wet" finally? Question is that during the water influx, the mined stope is like a unconfined aquifer? it is however in the middle, not the top.

Yes, my topic is an underground coal mine wtih exactly what you assumed. I have water samples from the stope before closure. But they are surprisingly not acidic(pH 7~8).  Probably they are mixed water and there may be buffering materials(back fillings) although we know pyrite does exist. I am trying to work out the process with PHAST. Many thanks for your detailed explanations and regards.

[dlparkhurst]

I don't think it is possible for PHAST to model a dry zone between two saturated zones. You can raise the water table to fill an unsaturated volume. Conceivably you could could use a flux boundary or well to add suitable mine water to the stope area and other flux, leaky, or constant head boundaries to add the background aquifer water. It is not possible to change boundary condition types, although a single face can have both flux and leaky boundaries.

Perhaps the filling process is not that important and you can start your simulation with a saturated system with coal-mine water in the stope.

[rapheul]

I also think Phast could not simulate such dry zones between aquifers at present. Perhaps yes, the filling process is not so important, at least if I want to deal with the problem using phast, I need to define it as a saturated zone. In fact, my idea is just to be clear with the flow field changes caused by closure at beginning and then could I apply it for transport later.
Ok, now I would like to say thanks again for your time and all kind help. that was really helpful to my topic.
Best regards   

[rapheul]

Hi David,

After a time of trial of my model, there are still multi layer problems(flow only), I always got error and exit saying that failure to convergence on steady state flow initial conditions. Therefore I guess the wrong parts should be initial water heads or rivers.

For differernt aquifers, Head_IC and -prism are used to define various initial heads, I used the same -prism which exactly media aquifers used to write initial water heads, and for aquitards, no initial heads assigned. There was always a default head for the whole domain though, even I tried with .dat files to write the Keyword block without it. Now I am not clear what would be the problem for model run. cause 1st kind BC is so far the only intial condition for the model. I did not assign wells bc or any other bc yet. I also tried to delete the rivers but it still saying the same. What is more strange, if you deactivate some layers, run only either a single aquifer, it is working. Followings are parts of the trans.dat file.

HEAD_IC
   -domain
      -head                     42
HEAD_IC
   -prism
      -description Aquifer Q
      -top       CONSTANT   GRID 50
      -bottom    CONSTANT   GRID 0
      -perimeter SHAPE      MAP C\Whole.1.shp
      -head                     42
HEAD_IC
   -prism
      -description Aquifer S
      -top       CONSTANT   GRID -20
      -bottom    CONSTANT   GRID -100
      -perimeter SHAPE      MAP C\Whole.1.shp
      -head                     -15

 
For river definitions, it is very strange when I extend the thickness of the unconfined aquifer to more than a certain value, it will work(single layer model). then if I shorten the Z range, the model is running with errors. But my first aquifer is not that thick although the whole area is a little bit large with a width of 4-5km. It seems to have something related with the Z range as it runs well when the range is large enough.

Hope this not bothering you too much as I have ran it with trial and error too many times and the quick way maybe is still asking someone. Thanks in advance.

[dlparkhurst]

If you are running a steady-flow simulation, the initial head distribution should not be critical. I would leave just the default definition of 42 (or perhaps some other value) everywhere; remove the others.

Not quite clear what you have for boundary conditions and rivers, but I would add a single fixed head boundary cell (head 42 for example), and remove all other boundary conditions and rivers. That should run.

You can then begin adding the other boundaries and rivers sequentially and try to determine the problem.

[rapheul]

Thanks for your reply, I did as you told me to remove all the others, but it still says errors. maybe you need the followings text for checking my definitions.
The two specified bc are for the confined aquifer, means a hydraulic gradient for driving the flow. I simplified the points coordinates to be space saving here.
And two rivers for unconfined aquifer. in between is an aquitard.

SPECIFIED_HEAD_BC
   -prism
      -top       CONSTANT   GRID -20
      -bottom    CONSTANT   GRID -100
      -perimeter POINTS     GRID
         X Y
         end_points
      -head
         0 seconds   -8
SPECIFIED_HEAD_BC
   -prism
      -top       CONSTANT   GRID -20
      -bottom    CONSTANT   GRID -100
      -perimeter POINTS     GRID
         X Y
         end_points
      -head
         0 seconds   -18
RIVER 1
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20497175.1698501 3797388.7901807
      -head
         0 days     48
      -width                       10
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20505011.4303025 3797680.69857402
   -point 20509318.333429 3796482.76059105
   -point 20509460.4980225 3792952.04624298
      -head
         0 days     48
      -width                       10
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
RIVER 2
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20516753.2428071 3804644.71164472
      -head
         0 days     38
      -width                       10
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20519211.9997946 3798674.29051045
      -head
         0 days     38
      -width                       10
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1

Not even once did it run well, each time no matter what I change or remove, it is still saying errors.

However, if I define the grid as only one single unconfined aquifer, sometimes it may run well and the result is also what I am expected. please see below. the differences between the two are only the thickness of the aquifer. but the first one which is thicker runs without errors and the second one fails. I am really confused about this. and the initial heads do infect the model run. e.g. the first one only 45 runs, other values do not.
1:
GRID
   -uniform X 20496000 20521000 51
   -uniform Y 3792000 3807000 31
   -uniform Z -100 50 11
   -snap X 0.001
   -snap Y 0.001
   -snap Z 0.001
   -chemistry_dimensions XYZ
   -print_orientation XY
   -grid_origin   0  0  0
   -grid_angle    0
MEDIA
   -domain
      -active                   0
      -Kx                       0.0001
      -Ky                       0.0001
      -Kz                       1e-005
      -porosity                 0.2
      -specific_storage         0.0001
      -long_dispersivity        1
      -horizontal_dispersivity  1
      -vertical_dispersivity    1
      -tortuosity               1
   -prism
      -description Aquifer Q
      -top       CONSTANT   GRID 50
      -bottom    CONSTANT   GRID -100
      -perimeter SHAPE      MAP C\PZWhole.1.shp
      -active                   1
      -Kx                       5
      -Ky                       5
      -Kz                       0.5
      -porosity                 0.2
      -specific_storage         0.01
HEAD_IC
   -domain
      -head                     45
CHEMISTRY_IC
   -domain
      -solution            1
RIVER 1
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20510583.7457693 3793574.83622841
      -head
         0 days     48
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20519132.4682666 3798827.62212079
      -head
         0 days     48
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
RIVER 2
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20502340.8760918 3804896.34474171
      -head
         0 days     40
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20500286.0883328 3800778.02592431
   -point 20498037.069831 3800635.40916317
      -head
         0 days     40
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1

2:
GRID
   -uniform X 20496000 20521000 51
   -uniform Y 3792000 3807000 31
   -uniform Z -100 50 11
   -snap X 0.001
   -snap Y 0.001
   -snap Z 0.001
   -chemistry_dimensions XYZ
   -print_orientation XY
   -grid_origin   0  0  0
   -grid_angle    0
MEDIA
   -domain
      -active                   0
      -Kx                       0.0001
      -Ky                       0.0001
      -Kz                       1e-005
      -porosity                 0.2
      -specific_storage         0.0001
      -long_dispersivity        1
      -horizontal_dispersivity  1
      -vertical_dispersivity    1
      -tortuosity               1
   -prism
      -description Aquifer Q
      -top       CONSTANT   GRID 50
      -bottom    CONSTANT   GRID -50
      -perimeter SHAPE      MAP C\PZWhole.1.shp
      -active                   1
      -Kx                       5
      -Ky                       5
      -Kz                       0.5
      -porosity                 0.2
      -specific_storage         0.01
HEAD_IC
   -domain
      -head                     45
CHEMISTRY_IC
   -domain
      -solution            1
RIVER 1
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20510583.7457693 3793574.83622841
      -head
         0 days     48
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20519132.4682666 3798827.62212079
      -head
         0 days     48
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
RIVER 2
   -xy_coordinate_system GRID
   -z_coordinate_system  GRID
   -point 20502340.8760918 3804896.34474171
      -head
         0 days     40
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1
   -point 20500286.0883328 3800778.02592431
   -point 20498037.069831 3800635.40916317
      -head
         0 days     40
      -width                       5
      -bed_hydraulic_conductivity  5
      -bed_thickness               1
      -depth                       1

[dlparkhurst]

Attach one .trans.dat file that can be read by Phast4Windows, and I will look at it. Please make it as simple as possible.

[rapheul]

OK, please find the attached file. Actually it is still a trying example for my topic as there are more aquifers contained in this large area. Unluckily it is not working but if I define with this grid only one single aquifer, sometimes it runs well even with these two rivers included. Anyway, thanks for your kind help~

[dlparkhurst]

I really don't think that  PHAST is capable of modeling dry zones beneath saturated zones.

However, I was able to run to steady state by running in the transient mode, using the following time sequence.

TIME_CONTROL
   -time_step
      0 0.1
      3.6 1
      36 10
      360 100
      3600 1000
      36000 10000
      360000 100000
   -time_change
      3650000
   -start_time 0

Note, I did not have the shp file, so I simply used the entire domain for the perimeter.