PHAST: output of concentrations in WELL doesn't seem to be correct

[Leo]

Hello fellow PHREEQC/PHAST users,

I am currently using PHAST to simulate processes of Uranium dissolution/precipitation. I added a well with a very small pumping rate to simulate an observation well.

The output of *.wel.xyz.tsv shows very high concentrations of U, compared to those visible in the model viewer visualization as well as the data written to *.chem.xyz.tsv (and *.h5) for the nodes around the well.

What can be the problem here or am I missing something? I had a look at the *.wel.txt which shows basically the same concentrations / mass fluxes as the .*wel.xyz.tsv file. I don't think these can be representative concentrations in the 'raw water' of the well.

Thanks in advance.

[dlparkhurst]

I don't know of a problem with the well calculation. The concentration should be the flow-averaged concentration, probably in mol/kgw.

[Leo]

I don't know of a problem with the well calculation. The concentration should be the flow-averaged concentration, probably in mol/kgw.

Thanks for the reply, I really appreciate it.

I understand the calculation of the flow-weighted concentrations in the well. But the mass flow rates - for each node in z-direction - are way too high. While the 'plume' does not arrive at the well location, the concentrations of the well water are already rising which seems impossible. I would like to show you some screenshots, but I am currently not in the office. Any other ideas?

[dlparkhurst]

I will need to see the input files. Hopefully as simple as possible.

Perhaps we should go offline if you do not want to post the files here.

[Leo]

I will set up a simple model to verify the problem and upload the input files tomorrow.

[Leo]

I managed to reproduce the problem with a simple model.

The 'pst.trans.dat'-file was exported from Phast4Windows. The database used is 'wateq4f.dat'.

If you compare the output of 'pst.wel.xyz.tsv' with the concentrations in the nodes where the well is placed (x=95, y=25) or with the model viewer visualization, you will notice an increasing U-concentration in the well water which should not be possible.

I also wanted to attach an Excel-file with the results but it's bigger than 128 kb, so I uploaded it to:
http://s000.tinyupload.com/?file_id=74435710624460391722

Thanks for your help.

I will need to see the input files. Hopefully as simple as possible.

Perhaps we should go offline if you do not want to post the files here.

[dlparkhurst]

I think the issue is that the wel.xyz.tsv file is written after the transport part of the operator splitting sequence, whereas the .h5 and .chem.xyz.tsv files are written after the chemistry part of the  sequence.  In your example, the U after transport is relatively large, but most of the U precipitates during the chemistry calculation, so the Model Viewer concentration is small. You can look at d_UO2(a), which is the amount of UO2 dissolved (positive) or precipitated (negative) during the chemistry calculation. If you plot it in Model Viewer, I think the d_UO2(a) is about the U concentration given in the wel.xyz.tsv file.

The difference is related to the operator splitting error, which is proportional to the time step. I think if you decrease the time step, the differences between the wel.xyz.tsv and the .h5 file will decrease.

[Leo]

That seems like a good explanation. I wanted to use the data from the '*.wel.xyz.tsv'-file to compare tracer/U breakthrough curves. In this case it might be better to look at the nodal output of '*.chem.xyz.tsv'. What I liked about the well-file is the composition of a weighted 'mixed' water.

Would it be possible to implement a feature like this in a future phast-version, maybe 'just' by calculating the well-chemistry after the chemistry step?

[dlparkhurst]

Ken Kipp was the one who knew the transport calculation and the well calculation. I am not likely to try to rewrite and debug the sequence of calculation and printing.

My suggestion would be to use a smaller time step after you are happy with the calculation. If you have a multiprocessor computer (or a cluster), the MPI version of the code should help with the run times. You can compare a sequence of smaller time steps to see how the breakthrough changes.

[Leo]

Using smaller time steps as you suggested, almost eliminates the problem. Also I could see, that UO2(a) from the *.chem.xyz.tsv fits the concentrations in *.wel.xyz.tsv as you mentioned.

Still, I will probably prefer using the nodal values for comparing the tracer breakthrough curves for my problem, because of the fine spatial discretization of the real model and therefore very long computation times. I do not have access to more powerful computers right now and using Phast4Windows already uses all of my four CPU cores at 100%. For the future, I will definitely keep all of this in mind.

Anyway, I want to thank you a lot for all of your help, not only for almost instantanously answering the questions on this forum, but for putting such effort in these programs, that really help us understand and model hydrogeochemical problems.

[dlparkhurst]

OK. There are some formulas in the appendix of the PHAST manual for numerical dispersion due to grid and time steps. There is not a good formula for operator splitting error, but you can estimate that with some experimentation with the time step.