PhreeqcUsers Discussion Forum

Processes => Reactive transport modelling => Topic started by: eng.youzerseef on February 07, 2020, 05:09:35 PM

Title: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 07, 2020, 05:09:35 PM

I want to model a column experiment where a high pH fluid is pumped from one side through a sandstone that contains Quartz, K-feldspar and Kaolinite. the fluid is collected from the other side from concentration analysis.
my aim is to model the rate of dissolution and precipitation of the minerals by kinetic datablock, then use the transport datablock to study the reaction in the whole column (I assume the column is one cell since it is filled with only sandstone ). because of the precipitation and dissolution process, the surface area, number of moles and porosity of each mineral will change and hence the total column porosity.  Hence, I want to incorporate variable surface area in the rate equation of each mineral such that it changes each step by the change in the porosity value, something like the below equation:

A0_new = A0 * (mineral_initial_Porosity / mineral_Porosity_2)^2/3   


-Do I have to use PhreeqcRM or I could do it in the regular Phreeqc ? if phreeqc is capable of doing this, how could I implement variable surface area and variable porosity ?


your corporation is highly appreciated...
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 10, 2020, 04:27:46 PM
If you want the change in porosity to affect advection, you need to use IPhreeeqc or PhreeqcRM. The only way to change advective velocity is to define a new TRANSPORT data block with a new time step, which would be done under program control.

Porosity can affect diffusive transport, and there are Basic functions CHANGE_POR and GET_POR that allow porosity to evolve under diffusion. You might be able to use these functions to store and retrieve porosity even if you are not simulating diffusion only.

If you only want to use porosity to adjust rates, then you should be able to do it in the RATES definition.  You would have to decide how you want to set it up, but Basic function EQUI and KIN could give the current moles of minerals present. EQUI_DELTA and KIN_DELTA give the change in moles from the current simulation.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: MichaelZ20 on February 10, 2020, 05:29:23 PM
See also
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 11, 2020, 03:18:12 PM

thanks a lot for your reply...

in the example below...there is only calcite in the reaction and hence the initial porosity is actually calcite porosity. my question is...since I have different mineral in the sandstone should I divide the initial porosity which was calculated from experiment to three minerals porosity...such that

Ini_Porosity= Quartz_porosity1 + Kfeldspar_Porosity1 + Kaolinite_porosity1

and then the new porosity will be calculated for each mineral (based on the new volume which can be calculated from the volume of the dissolved or precipitated moles) and then...

New_Porosity=Quartz_porosity2 + Kfeldspar_Porosity2 + Koalinite_porosity2

of this approach is fine I am facing another problem which is the reliability of the PUT & GET functions since it is mention in the manual that those functions are unreliable in parallel processes.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 11, 2020, 03:33:20 PM
I am revising this post. GET and PUT have limitations for both OpenMP and MPI parallelization.

PUT and GET should work for passing information within a cell at a given time step; for example from one RATES definition and another for the the same cell during the same time step. The values are not guaranteed between cells or between time steps. The PUT values are stored on the processor (or thread) where the PUT is executed, but cells are distributed among all the processors (or threads) (and may be transferred to balance the load on each processor or thread). So a GET statement could be executed on a different processor (or thread) from the one the PUT statement was executed on.

Title: Re: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 12, 2020, 03:34:04 PM
Thanks Dr. Parkhurst

what about my point of dividing the total porosity among the minerals precipitants in the sandstone sample. do you think it is valid?
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 12, 2020, 03:53:06 PM
What you really want is reactive surface area, but even then the rate of reaction can depend on the number and type of defects on the mineral surface. Most phreeqc.dat rate equations simply use the (M/M0)^n, where M is current moles, M0 is initial moles, and n is an exponent. Seems like you are doing almost the same thing. With all the uncertainties, I don't think one approach clearly better.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 21, 2020, 05:10:36 PM

i am trying to create a loop in Phreeqc for the transport time steps but it seems like it is not working.

my rate equation is:

rate_Q = A0_Q * (QMass) *  ((m/m0)^0.67) * (10^-12.7) * (1 - SR("Quartz"))

-reactive surface area are in m2/g, that why i multiple by the mass (Qmass)

-I need to change the below variables (The first two variables in the rate equation) starting from the second run.
               #A0_Q= Parm(1) for the first run, and then A0_Q = Parm(1) * (Por1 / Por2)^2/3 
               #QMass=Qmass1 for the first run, then Qmass = Qmass1 - (Quartz_dissolve _moles_mass)
-I want the program to read the normal rate equation(with initial reactive surface area (parm(1)) and initial mass) for the first run and then based on the new variable the rate equation should be modified.

- I tried the below:

If SIM_TIME=0 then A0_Q = Parm(1) else A0_Q = parm(1) * (Por1 / Por2)^2/3 

IF SIM_TIME=0 then QMass=Qmass1 else QMass=Qmass1-(Quartz_dissolve _moles_mass).

- I have two problem with the results:
1- the program start from the second time step and not 0 (the results of time=0 are shown in the second time step and it is all zeros at time =0)
2- the program work fine for the first and second run, but then the values are all constant as the second run values.

any ideas what could be done?
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 21, 2020, 11:41:29 PM
You don't say how you are calculating time, but TOTAL_TIME is not equal to zero during the first time step. The rates are evaluated at a number of times over the first time_step interval, so TOTAL_TIME is probably zero only for one sub time step of the integration.

If M = M0, I don't see why A0_Q = parm(1) * (Por1 / Por2)^2/3  would not work for all time. Again, as the integration proceeds over a single time step, porosity would vary as estimated at different TOTAL_TIME (over the time step interval) to achieve an overall rate for the entire time step. So, the change in porosity (Por1/Por2) is not constant, but evaluated at each of multiple sub time steps. In general, continuous functions should be used for the rate equations, including continuity within a given time step interval.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 22, 2020, 04:18:51 PM

-I have the total time for the conducted experiment which is 5590 hours. in my modelling i used 26 hour time step with 215 shifts. i did not include any time steps in the kinetics data black.

-what do you mean by continuous equation for the rate.

- I have attached my phreeqc file.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 22, 2020, 07:38:10 PM
TRANSPORT uses the time step defined in the transport datablock, regardless of the KINETICS definition.

You should add some END statements to your file. PHREEQC does a default reaction calculation if a SOLUTION is defined within the same block demarcated by END statements. In your file SOLUTION 0 is reacted with EQUILIBRIUM_PHASES 0 and KINETICS 1. The result of the reaction is not SAVEd, so the SOLUTION 0 and EQUILIBRIUM_PHASES 0 are not modified. KINETICS 1 is updated, so the moles of kinetic reactant is changed before the TRANSPORT calculation. I would probably add an END after each keyword data block.

The rate is probably continuous if you do not have IF statements. Your original description talked about doing something different for the first step than other steps, which could result in a discontinuity in the rate equation.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: eng.youzerseef on February 24, 2020, 04:12:59 PM
Thank you so much prof. Parkhurst.

I tried to add END in my input file, but after running...the program only doing initial solution and transport.

If I want the kinetics reaction to perform I have to remove all END from solution 0, solution 1-10 and equilibrium phase (as my attached file).

I even tried to run example 15 (In the software it runs without END for solution 0 and solution 1-10. I found an updated version of the example online which have END after solution 0 and solution 1-10), but in both cases the same thing is happening (No kinetic reaction, only transport).

I guess I still do not fully understand how the program works and when should I save the solution for the next calculation.

anyway, My reaction should be like this:

1-background fluid with sandstone
2-pulse fluid is pumped, kinetic reaction take place during the transport of the pulse fluid.

i really appreciate all your support.
Title: Re: Column Exp. & Reactive Transport Modelling
Post by: dlparkhurst on February 24, 2020, 06:03:29 PM
If you do not understand the calculations, then you should look at the output file. Below is a synopsis of the output generated by your input file. The sequence of initial, reaction, transport, and other calculations is described in the manual in the section Types and Sequence of Calculations.

SAVE applies only to the "batch-reaction calculations". To simplify, I suggest you use an END after each keyword data block. If you put END after each keyword data block in your file, you will not get a batch-reaction calculation.

If you want to react one of the solutions with reactants and save results to initialize the TRANSPORT calculation, I think you should define the batch reaction explicitly. In your case, to recreate your output, you would add the following after KINETICS ... END:

USE solution 0
USE equilibrium_phases 1
USE kinetics 1

This is the reaction defined by default in your original input file. If you want to save the solution created by this reaction as the new solution 0 for TRANSPORT, you would add

SAVE solution 0 as part of the explicit reaction definition.

If you want to save the new moles of equilbrium phases, add SAVE equilibrium_phases n as part of the explicit reaction calculation, where n is probably 0 or 1.

Unlike solution or equilibrium phases, the reaction (whether the default defined in your original file or the explicit reaction defined above) will automatically update the moles of the kinetic reaction of KINETICS 1, so that the updated moles will be used in cell 1 in the TRANSPORT calculation. If you add the ENDs and do not define the explicit reaction, the moles of KINETIC reaction used in TRANSPORT will be as originally defined in KINETICS.

So, yes, the different input files have a kinetic batch reaction or not. Your choice.

Beginning of initial solution calculations.

Initial solution 0.   Pulse solution with YCL

Initial solution 1.   Background solution initially filling column

Beginning of batch-reaction calculations.

Reaction step 1.

WARNING: Element Al is contained in K-Feldspar (which has 0.0 mass),   
but is not in solution or other phases.
Using solution 0.   Pulse solution with YCL
Using pure phase assemblage 1.   
Using kinetics 1.   

Reaction step 2.

Reaction step 215.

Reading input data for simulation 2.
Beginning of transport calculations.