Conceptual Models > Incorporation PHREEQC in programming languages


(1/3) > >>

Dear all,

I am working on coupling PhreeqcRM with a transport model written by Fortran. In advection_f90.F90, there are several sentences attached below. I know the model domain is symmetric in the given example, but the concentration of each cell is not. So, why can we treat this as two equivalent rows by symmetry?

--- Code: ---  ! Demonstation of mapping, two equivalent rows by symmetry
  do i = 1, nxyz/2
     grid2chem(i) = i - 1
     grid2chem(i+nxyz/2) = i - 1
--- End code ---

In addition, I don't understand the difference between workers, initial_phreeqc and utility. As I understand, I first transfer the initial condition from PhreeqcRM to the transport model and hence I should set initial_phreeqc as 1 and others as 0. During the transport calculation, I should set workers as 1 and others as 0. Am I correct?

Any explanation would be appreciated.

Best regards,

Just for example, the grid for HST3D is a rectangular point-distributed grid. There are at least two nodes in each direction, so consider a 1-dimensional calculation that has two nodes in Y and Z and 10 nodes in X. There are basically four columns of cells in the X direction, and if initial and boundary conditions are the same for all four, then the flow and chemistry will evolve identically. You could calculate chemistry for all 40 cells or you could calculate chemistry for one column and copy it to the other three columns. PHAST, in fact, allows this option to speed up the chemistry calculation for 1D problems. If that is not clear to you, then let's skip it and you can consider only non-symmetric flow and chemistry.

The way PhreeqcRM is designed to work is to define enough solutions and other chemical reactants (like EQUILIBRIUM_PHASES) in the Initial IPhreeqc instance to be able to define all the initial conditions for the model. The user numbers for the solutions and reactants are used in RM_InitialPhreeqc2Module to populate all the model cells, which are in the worker instances. The Initial IPhreeqc instance is also used to define boundary-conditions solutions; the concentrations of the components of these solutions can be retrieved by RM_InitialPhreeqc2Concentrations for boundary conditions in your transport model. So the initial iphreeqc instance is a work space where you can run Phreeqc calculations to obtain initial and boundary conditions for the model. You can do these calculations by setting only the initial iphreeqc flag to 1 in RM_RunFile and RM_RunString.

However, if you define SELECTED_OUTPUT and USER_PUNCH, you want make sure this is defined to the workers so that their output is available as the reactive-transport calculations proceed. Also, thermodynamic information, like SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, PHASES, SURFACE_SPECIES and others, should be run for the workers, initial iphreeqc, and utility so that all of the IPhreeqc instances have the same thermodynamic model.

After initial conditions are established in the workers for all nxyz cells, then you may want to use one RM_RunCells with a zero time step to establish equilibrium between the solutions and other reactants of the model. At that point the chemistry is initialized and you can use RM_GetConcentrations to obtain the concentrations of all the components for all of the cells of your model for your first transport time step. At this point, you probably will not reset any initial conditions and you will not use RM_InitialPhreeqc2Module again. You may need to use RM_InitialPhreeqc2Concentrations to obtain boundary conditions for subsequent changes in boundary conditions.

The utility IPhreeqc instance is seldom used and you can likely ignore it. In PHAST, the utility is used for wells. The water that is extracted from a well is a mixture of concentrations from multiple model cells. The number of moles of each component are known for the extracted water, but the pH and, say, saturation indices for minerals for the water are not known. The extracted well-water composition is sent to the utility IPhreeqc instance with RM_Concentrations2Utility and is then processed to calculate pH and other properties.

Dear Parkhurst,

Thank you very much for your reply.

For now, I have separated advection.f90 from the whole project and ran it successfully by including RM_interface.f90. In this way, I can check each parameter by printing it on the screen. Inspired by this example, I figured out how to implement PhreeqcRM to my transport model.

However, I have another two questions. The first one is how to set the pressure in the unsaturated zone for PhreeqcRM. Since the transport model calculates the pressure head in both the saturated and unsaturated zones, how can I transfer this pressure head to the pressure used in PhreeqcRM? In addition, I don't understand the sequence of the components after reaction. Is there any rule to define this component sequence?

By the way, I may calculate pe values through this coupling model. So, I was wondering what extent can I trust the calculated pe values from PhreeqcRM?

Thanks for your help in advance.

Best regards,

If you are asking how calculated negative pressure heads affect geochemical reactions, I have no idea. You will have to look in the literature. My simple mind would set the pressures of unsaturated cells to 1 atm. I can say that the pressure dependence of equilibrium constants in the Phreeqc databases  are fit for 1 atm and above for those databases that have pressure dependence.

Look at the documentation for RM_FindComponents, RM_GetComponentCount, and RM_GetComponent. The order of the components is given by calling RM_GetComponent in sequence after RM_FindComponents has been called.

Search the forum for "redox equilibrium" for discussions of pe.

Dear Parkhurst,

Thanks for your reply.

I have looked at the documents you mentioned. I know I can get a list of components after calling RM_GetComponent. Howerver, I am still questioned about the order of components in that list, since I have to define the boundary condition and initial condition correspondingly for each component in the transport model eralier. If I don't know the order of conponents earilier, how do I set corresponding boundary condition and initial condition for each component in the transport model? Thanks a lot.

Best regards,


[0] Message Index

[#] Next page

Go to full version