Version 3.8.0-16920
Version 3.8.0: July 11, 2024
-----------------
July 11, 2024
-----------------
PHREEQC: Fixed a bug in the DUMP routines. Under some circumstances
erroneous output was dumped for a user number. In most cases, the
correct output was dumped following the erroneous output, which
caused the erroneous output to be ignored.
-----------------
May 18, 2024
-----------------
DATABASES:
sit.dat was updated to version 12a (Aug 22, 2023) from www.thermochimie-tdb.com.
Amm.dat, iso.dat, llnl.dat, minteq.dat, minteq.v4.dat, phreeqc.dat,
phreeqc_rates.dat, pitzer.dat. Tipping_Hurley.dat, and wateq4f.dat were
reformatted by using the lsp utility by David Kinniburgh from phreeplot.org.
-----------------
May 3, 2024
-----------------
PHREEQC: The -dw identifier of SOLUTION_SPECIES now has up to 7 items.
-dw Dw(25C) dw_T a a2 visc a3 a_v_dif
where,
Dw(25C)--Tracer diffusion coefficient for the species at 25 °C, m 2 /s.
dw_T--Temperature dependence for diffusion coefficient.
a--Debye-Hückel ion size.
a2--exponent.
Visc--Viscosity exponent.
a3--Ionic strength exponent.
A_v_dif--Exponent for (viscosity_0/viscosity).
The diffusion coefficient is calculated as follows:
Dw = Dw(25C) * exp(dw_T / T - dw_T / 298.15)
ka = DH_B * a2 * I0.5/ (1 + a3)
av = (viscos_0/viscos)a_v_diff
ff = av * exp(-a * DH_A * z * I0.5 / (1 + ka))
Dw = Dw * ff
Where T is temperature in Kelvin, DH_B is the Debye-Hückel B parameter,
I is ionic strength, viscos_0 is the viscosity of pure water at T, viscos is
the viscosity of the solution at T, DH_A is the Debye-Hückel A parameter,
and z is the charge on the species,the viscosity of the solution.
See Robinson and Stokes, 2002, Chpt 11 for examples.
The Dw and a_v_dif can be set in a USER_ program with
setdiff_c("name", Dw, a_v_dif), for example:
10 print setdiff_c("H+", 9.31e-9, 1).
The diffusion coefficient of H+ is handled differently with
Falkenhagen equations.
-----------------
May 3, 2024
-----------------
PHREEQC: The ionic strength correction is for electromigration calculations
(Appelo, 2017, CCR 101, 102). The correction is applied when the 6th parameter
option is set to true for -multi_D in TRANSPORT:
-multi_d true/false 1e-9 0.3 0.05 1.0 true/false # multicomponent diffusion
true/false, multicomponent diffusion is used,
default tracer diffusion coefficient (used in case -dw is not defined for a species),
porosity (por = 0.3),
limiting porosity (0.05) below which diffusion stops,
exponent n (1.0) used in calculating the effect of tortuosity on the
porewater diffusion coefficient Dp = Dw * por^n,
true/false: correct Dw for ionic strength (false by default).
-----------------
May 3, 2024
-----------------
Database: Added new database phreeqc_rates.dat. The database augments
phreeqc.dat with rate parameters from Palandri and Kharaka (2004),
Sverdrup, Oelkers, Lampa, Belyazid, Kurz, and Akselsson (2019) (only
Albite and quartz), and Hermanska, Voigt, Marieni, Declercq,
and Oelkers (2023). Parameters are defined in data blocks
RATE_PARAMETERS_PK, RATE_PARAMETERS_SVD, and RATE_PARAMETERS_HERMANSKA.
All minerals with rate parameters have been added in a PHASES
data block. Example RATES definitions using the different RATE_PARAMETERS_
parameters are provided for Albite and Quartz.
-----------------
April 27, 2024
-----------------
Databases: Added new keyword data block MEAN_GAMMAS. Each line
of the data block defines how to calculate the mean activity
coefficient for a salt with a series of pairs of
aqueous species and stoichiometric coefficient. Phreeqc.dat,
Amm.dat, pitzer.dat, and phreeqc_rates.dat have this data block.
MEAN_GAMMAS
MgCl2 Mg+2 1 Cl 2
A new Basic function MEANG will calculate mean activity coefficients
for salts listed in the MEAN_GAMMAS data block.
10 g_MgCl2 = MEANG("MgCl2")
-----------------
April 27, 2024
-----------------
PHREEQC: Added new keyword data blocks RATE_PARAMETERS_PK, RATE_PARAMETERS_SVD,
and RATE_PARAMETERS_HERMANSKA and Basic functions RATE_PK, RATE_SVD, and
RATE_HERMANSKA
RATE_PARAMETERS_PK
# Acid Neutral Base
# log K E n(H+) log K E log K E n(OH-)
# ======== ======== ======== ======== ======== ======== ======== ========
Quartz -30 0 0 -13.4 90.9 -30 0 0 # Table 4
# Acid Neutral P_CO2
# log K E n(H+) log K E log K E n(P_CO2) Table
# ======== ======== ======== ======== ======== ======== ======== ======== ========
calcite -0.3 14.4 1 -5.81 23.5 -3.48 35.4 1 33 # specify Table number for P_CO2^n(P_CO2)
# Acid and Fe+3 Neutral and O2 Base
# log K E n(H+) n(Fe+3) log K E n(O2) log K E n(OH-) Table
# ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ========
pyrite -7.52 56.9 -0.5 0.5 -4.55 56.9 0.5 -30 0 0 35 # specify Table number for Fe+3 and O2
Three rate equations from Palandri and Kharaka (2004) can be entered. Most minerals use
use the first form above with 8 parameters. Table 33 has a term for CO2 as in
the calcite example above; parameters from table 33 are identified with a 33 in the 9th
field following 8 parameters. Table 35 has additional terms and data from this table
is identified with 35 in field 11 following 10 rate parameters. The rates for the
the minerals listed in the data block can be calculated with the Basic function RATE_PK.
The calculated rate does not include factors for surface area or affinity.
10 rate = RATE_PK("Calcite")
RATE_PARAMETERS_SVD
# Table 4: E's Table 3: H+-reaction H2O-reaction CO2-reaction Organic_acids OH--reaction Table 5
# H+ H2O CO2 Org_acids OH- pkH nH yAl CAl xBC CBC pkH2O yAl CAl xBC CBC zSi CSi pkCO2 nCO2 pkOrg nOrg COrg pkOH- wOH- yAl CAl xBC CBC zSi CSi # Num Mineral Formula
# ====== ====== ====== ========= ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ===== ====== ====== ====== ====== ====== ====== ====== ====== ====== ====== ======= ======
Albite 3350 2500 1680 1200 3100 14.6 0.5 0.4 0.4 0.4 0.5 16.8 0.15 4 0.15 200 3 900 16.05 0.6 14.7 0.5 5 15.4 0.3 0.1 12 0.5 5 3 900 # 1.6 Albite NaAlSi3O8
Rate parameters from Sverdrup, Oelkers, Lampa, Belyazid, Kurz, and Akselsson (2019)
can be specified with the RATE_PARAMETERS_SVD data block. A total of 31 parameters
are entered for each mineral. The rates for minerals minerals listed in the data
block can be calculated with the Basic function RATE_SVD. The calculated rate does
not include factors for surface area or affinity.
10 rate = RATE_SVD("Albite")
RATE_PARAMETERS_HERMANSKA
# Acid mechanism Neutral mechanism Basic mechanism
# logk25 Aa Eaa n(H+) logk25 Ab Eab logk25 Ac Eac n(OH) # Formula
# ======== ========= ======== ======== ======== ========= ======== ======== ========= ======== ======== =========================================
# Amphiboles
Anthophyllite -12.4 5.70E-04 52 0.4 -13.7 5.00E-06 48 0 0 0 0
Rate parameters from Hermanska, Voigt, Marieni, Declercq, and Oelkers (2023) can
be specified with the RATE_PARAMETERS_HERMANSKA data block. A total of 11 parameters
are entered for each mineral. The rates for minerals listed in the data block can
be calculated with the Basic function RATE_HERMANSKA. The calculated rate does not
include factors for surface area or affinity.
10 rate = RATE_HERMANSKA("Anthophyllite")
-----------------
April 21, 2024
-----------------
PHREEQC: Added Basic functions GET$ and PUT$. They are are the same as
GET and PUT, except the first argument for PUT$ is a character string,
and GET$ returns a character string. You may use one or more indices as
needed to identify the value that is saved (PUT$) or retrieved (GET$).
PUT$("MgCl2", 1, 1, 1)
x$ = GET$(1, 1, 1)
-----------------
April 19, 2024
-----------------
DATABASE: Kinec.v2.dat is a new llnl.dat style database from the
CarbFix2 and GECO projects that is included in new distributions of
PHREEQC. This database contains the parameters for calculating mineral
dissolution rates for primary and secondary silicate minerals using the
equations and parameters reported by Hermanska et al. (2022, 2023)
and dissolution rates for other non-silicate mineral systems using the
equations and parameters reported by Oelkers and Addassi (2024, in
preparation).
-----------------
April 15, 2024
-----------------
PHREEQC: Fixed a memory error with iso.dat because it uses H3O+ instead of
H+. The SC variable was uninitialized in that situation.
DATABASES: Amm.dat, phreeqc.dat, and pitzer.dat were updated with
revisions to viscosity and specific conductance.
PhreeqcRM and IPhreeqc: Fixed bug with the temperature grid for llnl. Some
internal testing and list generators used the default temperature of 25C,
which caused an error if the temperature grid did not span 25C.
-----------------
March 25, 2024
-----------------
DATABASES phreeqc.dat, Amm.dat, and pitzer.dat: The calculation of the
specific conductance can now be done with a Debye-Hückel-Onsager equation
that has both the electrophoretic and the relaxation term. (The standard
phreeqc calculation uses a simple electrophoretic term only.) For
individual ions, the equation can be multiplied with the viscosity ratio of
the solvent and the solution, and the ion-size a in the Debye-Hückel term
kappa_a can be made a function of the apparent molar volume of the ion. The
options are described and used in the databases. The additions extend the
applicability of the DHO equation to concentrations in the molar range,
reducing AARD (average of the absolute relative deviations) for SC and
transference numbers to less than 1% in many cases. For high KHCO3
concentrations, the SCs indicate the presence of a KHCO3 complex that was
added to phreeqc.dat and Amm.dat. The AARD's are 0.18 % for NaCl, 0.48 %
for KCl, 0.51 % for MgCl2 and 0.89 % for CaCl2. More example files are
available at http://hydrochemistry.eu.
PHREEQC Bug-fix: Option -density c[alculate] in SOLUTION_SPREAD was
corrected to give the iterated density of the solutions.
PHREEQC: A new option has been added. The viscosity of the EDL
layer on SURFACE(s) can now be calculated and will then be used to
modify the diffusion coefficients. It is set by adding c(alculate)
after viscosity, for example, "-donnan 1e-8 viscosity calc".
PHREEQC Bug-fix: Viscosity of the EDL layer on SURFACE(s), defined with, for
example, "-donnan 1e-8 viscosity 3", was omitted in Version 3.4.2. It is
now re-introduced in the calculations.
PHREEQC Bug-fix: Basic now returns the contributions to the specific conductance
(t_sc("H+")) and the viscosity (f_visc("H+")) only when the species is present
in the solution. In previous versions a dummy value was returned when the
species was predefined, but absent in the actual solution calculation.
PHREEQC Bug-fix: Limits for fugacity coefficients were set to be 0.01 < phi < 85 in
Peng-Robinson calculations. The limits were removed in version 3.7 (when calculating
H2S(g) solubilities). However, without the limits, all water turned into H2O(g) in some
cases and calculations failed.
-----------------
November 15, 2023
-----------------
PHREEQC programs: Fixed a couple malloc checks, some compiler warnings, and removed
some deprecated calls to strcpy and strcat.
-----------------
November 5, 2023
-----------------
PHREEQC programs: Automatic testing was expanded to include MPI and additional compilers.
-----------------
November 1, 2023
-----------------
PHREEQC: Logical statement in k_temp was modified to work with Intel optimization.
The statement at the beginning of the routine was not handled correctly when some
values were NaN.
-----------------
August 29, 2023
-----------------
PhreeqcRM: Fixed bug in memory allocation for selected output. One array accumulated lines
indefinitely, leading to ever increasing memory use. Memory use should now be relatively
constant once all selected output has been defined and used.
-----------------
June 1, 2023
-----------------
PhreeqcRM: Finalizing a Python version of PhreeqcRM that includes the BMI capabilities.
Methods are documented in Python style and two test cases are available, one
of which uses every Python method that is available.
-----------------
May 22, 2023
-----------------
PhreeqcRM: Revised all F90 methods that return arrays to use allocatable arrays,
so that, getter arrays are automatically dimensioned to the correct sizes
-----------------
May 22, 2023
-----------------
PHREEQC: (See https://hydrochemistry.eu/ph3/release.html for html version of changes.)
Added Basic function f_visc("H+") that returns the fractional contribution of a species to
viscosity of the solution when parameters are defined for the species with -viscosity.
Actually, it gives the contribution of the species to the B and D terms in the Jones-Dole
eqution, assuming that the A term is small. The fractional contribution can be negative, for
example f_visc("K+") is usually less than zero.
Bug-fix: When -Vm parameters of SOLUTION_SPECIES were read after -viscosity parameters, the
first viscosity parameter was set to 0.
Defined -analytical_expression and -gamma for Na2SO4, K2SO4 and MgSO4 and Mg(SO4)2-2 species in
phreeqc.dat and Amm.dat, fitting the activities from pitzer.dat from 0-200 °C, and the solubilities of
mirabilite/thenardite (Na2SO4), arcanite (K2SO4), and epsomite, hexahydrite, kieserite (MgSO4
and new species Mg(SO4)2-2). The parameters for calculating the apparent volume (-Vm) and the
diffusion coefficients (-Dw) of the species were adapted using measured data of density and
conductance (SC). Example files are available at http://hydrochemistry.eu
Removed the NaCO3- species in PHREEQC.dat since it is not necessary for the calculation of
the specific conductance (SC) and its origin is unknown.
Defined parameters in the -analytical_expression, -gamma, -dw, -Vm and -viscosity for
the NaHCO3 species in phreeqc.dat and Amm.dat, using the data in Appelo, 2015, Appl. Geochem.
55, 62-71. (These data were used for defining interaction parameters in
pitzer.dat.) The parameters for the apparent volume (-Vm), the diffusion
coefficient (-Dw) and the viscosity of CO3-2 and HCO3- were adapted using measured
data of density, conductance and viscosity of binary solutions.
The viscosity of the solution at P, T is now calculated and printed in the output file, and can
be retrieved in Basic programs with the function viscos (in previous versions, viscos returned
the viscosity of pure water at P, T).
The calculation uses a modified Jones-Dole equation which sums the contributions of individual
solutes:
eta / eta0 = 1 + A sqrt(0.5 sum(zi*mi)) + sum fan (Bi*mi + Di*mi*ni),
where eta is the viscosity of the solution (mPa s), eta0 is viscosity of pure water at the
temperature and pressure of the solution, mi is the molality of species i, made dimensionless
by dividing by 1 molal, and zi is the absolute charge number. A is derived from Debye-Hückel
theory, and fan, B, D and n are coefficients that incorporate volume, ionic strength and
temperature effects.
The coefficients are:
B = b0 + b1 exp(-b2 tC)
where b0, b1, and b2 are coefficients, and tC is the temperature in °C. The temperature is
limited to 200 °C.
fan = (2 - tan * Van / VCl-)
for anions, with tan a coefficient and Van the P, T and I dependent, apparent volume of the
anion relative to the one of Cl-, which is used as reference species. For cations, fan = 1
and tan need not be defined.
D = d1 exp(-d2 tC)
where d1 and d2 are coefficients.
n = ((1 + fI)^d3 + ((zi^2 + zi) / 2 * mi)^d3 / (2 + fI)
where fI averages ionic strength effects and d3 is a coefficient.
The coefficients are fitted on measured viscosities of binary solutions and entered
with item -viscosity under keyword SOLUTION_SPECIES, for example for H+:
SOLUTION_SPECIES
H+ = H+
-viscosity 9.35e-2 -8.31e-2 2.487e-2 4.49e-4 2.01e-2 1.570 0
# b0 b1 b2 d1 d2 d3 tan
When the solute concentrations are seawater-like or higher, the viscosity is different
from pure water (see figure at). To obtain a valid model for natural waters with phreeqc.dat,
the complexes of SO42- with the major cations were redefined, as noted above.
The A parameter in the Jones-Dole equation needs temperature dependent diffusion coefficients
of the species, and therefore the parameters for calculating the I and T dependency of the
diffusion coefficients (-dw parameters of SOLUTION_SPECIES) were refitted for SO42- and CO32-
species. Example files are available at http://hydrochemistry.eu.
Implicit calculations with option -fix_current will now account for changing concentrations in
the boundary solutions of the column.
Activated the print of statements defined in USER_PRINT when the initial EXCHANGE, SURFACE and
GAS_PHASE are calculated.
Changed the dw_t parameter for CO3-2 to 30 (was 0) and for HCO3- to -150 (was 0) to better fit
McCleskey's data
Bug fix: removed the factor (TK / 298.15) from the calculation of the temperature dependence of
the diffusion coefficient. For an example, see the calculation of Dw(TK) of H+ in the next
paragraph.
Bug fixes in printing/punching of diffusion coefficients with diff_c and setdiff_c: the numbers
are now corrected for I and T effects when the appropriate factors are defined in keyword
SOLUTION_SPECIES, item -dw. For example:
H+ = H+
-gamma 9.0 0
-dw 9.31e-9 1000 0.46 1e-10 # The dw parameters are defined in Appelo, 2017, CCR 101, 102-113.
It will set Dw(TK) = 9.31e-9 * exp(1000 / TK - 1000 / 298.15) * viscos_0_25 / viscos_0_tc
and Dw(I) = Dw(TK) * exp(-0.46 * DH_A * |zi| * I 0.5 / (1 + DH_B * I 0.5 * 1e-10 / (1 + I 0.75))),
where viscos_0_25 is the viscosity of pure water at 25 °C, viscos_0_tc is the viscosity of pure
water at the temperature of the solution. DH_A and DH_B are Debye-Hückel parameters,
retrievable with PHREEQC Basic.
The temperature correction is always applied in multicomponent, diffusive transport and for
calculating the viscosity.
The ionic strength correction is for electromigration calculations (Appelo, 2017, CCR 101, 102).
The correction is applied when the option is set true in TRANSPORT, item -multi_D:
-multi_d true 1e-9 0.3 0.05 1.0 true # multicomponent diffusion
# true/false, default tracer diffusion coefficient (Dw = 1e-9 m2/s) in water at 25 °C (used in
case -dw is not defined for a species), porosity (por = 0.3), limiting porosity (0.05) below
which diffusion stops, exponent n (1.0) used in calculating the porewater diffusion coefficient
Dp = Dw * por^n, true/false: correct Dw for ionic strength (false by default).
-----------------
May 19, 2023
-----------------
PhreeqcRM:
Changed documentation of GetDensity and related functions to GetDensityCalculated.
(GetDensity still exists for backward compatibility.)
Changed documentation of SetDensity and related functions to SetDensityUser.
(SetDensity still exists for backward compatibility.)
Density is used to convert user-model concentrations to module solution definitions only if the
units of the user-model concentrations are specified to be parts per million. The density specified by
SetDensityUser is used by SetConcentrations to convert from per kg of solution to
per L of solution. For GetConcentrations, two options are available to convert from module solutions
to user-model concentrations, depending on the value used for the method SetUseSolutionDensityVolume:
(1) the module-calculated density is used to convert from the calculated volume of solution
to the mass (kg) of solution, or (2) the user-specified value of density is used to make the conversion.
Again, density is only used if the user-model concentration units are ppm.
The change in method names is intended to emphasize the difference between the user-specified densities
and the module-calculated densities.
Changed documentation of GetSaturation and related functions to GetSaturationCalculated.
(GetSaturation still exists for backward compatibility.)
Changed documentation of SetSaturation and related functions to SetSaturationUser.
(SetSaturation still exists for backward compatibility.)
The values specified by SetSaturationUser are used to convert user-model concentrations to module solution definitions.
For SetConcentrations, the volume of solution is calculated to be the user-specified saturation * porosity *
representative volume. For GetConcentrations, two options are available to determine the solution volume, depending
on the value specified for SetUseSolutionDensityVolume: (1) the solution volume is calculated by the reaction module
and used to convert to user-model concentrations, or (2) the solution volume is again calculated by
user-specified saturation * porosity * representative volume, and those values are used to convert to user-model
concentrations. In either case, the values returned by GetSaturationCalculated are the calculated solution volume divided
by (porosity * representative volume).
The change in method names is intended to emphasize the difference between the user-specified saturations and
and the module-calculated saturations.
-----------------
April 16, 2023
-----------------
PhreeqcRM: Added new methods to simplify getting and setting component and
aqueous species concentrations.
New methods:
GetIthConcentration(int i, std::vector<double>& c)--Gets the ith component concentration as
of the last RunCells calculation. Total number of components is retrieved with GetComponentCount.
GetIthSpeciesConcentration(int i, std::vector<double>& c)--Gets the ith aqueous species concentration as
of the last RunCells calculation. Total number of aqueous species is retrieved with GetSpeciesCount.
This method is for use with multicomponent diffusion, and SetSpeciesSaveOn must be set to true.
SetIthConcentration(int i, std::vector<double>& c)--Sets the ith component concentration; done after
transport calculations and before RunCells calculation. Total number of components is retrieved
with GetComponentCount. SetIthConcentration must be run for every component before RunCells is called.
SetIthConcentration(int i, std::vector<double>& c)--Sets the ith aqueous species concentration; done after
transport calculations and before RunCells calculation. Total number of aqueous species is
retrieved with GetSpeciesCount. This method is for use with multicomponent diffusion,
and SetSpeciesSaveOn must be set to true. SetIthSpeciesConcentration must be run for every aqueous
species before RunCells is called.
Fortran versions are
RM_GetIthConcentration(id, i, c)
RM_GetIthSpeciesConcentration(id, i, c)
RM_SetIthConcentration(id, i, c)
RM_SetIthSpeciesConcentration(id, i, c)
-----------------
April 14, 2023
-----------------
PhreeqcRM: Added new methods to simplify setting initial conditions.
New initial conditions methods:
InitialEquilibriumPhases2Module(equilibrium_phases);
InitialExchanges2Module(exchanges);
InitialGasPhases2Module(gas_phases);
InitialKinetics2Module(kinetics);
InitialSolutions2Module(solutions);
InitialSolidSolutions2Module(solid_solutions);
InitialSurfaces2Module(surfaces);
These methods are an alternative to InitialPhreeqc2Module, which used a 7 x nxyz array to
set all initial conditions in one method call. The new methods set only one reactant at a
time, and all methods use a single array of index numbers (referring to definitions in
the InitialPhreeqc instance) of length nxyz (the number of user grid cells). The methods
copy definitions from the InitialPhreeqc instance to define initial conditions in the
model cells.
Fortran implementation is as follows:
RM_InitialEquilibriumPhases2Module(id, equilibrium_phases);
RM_InitialExchanges2Module(id, exchanges);
RM_InitialGasPhases2Module(id, gas_phases);
RM_InitialKinetics2Module(id, kinetics);
RM_InitialSolutions2Module(id, solutions);
RM_InitialSolidSolutions2Module(id, solid_solutions);
RM_InitialSurfaces2Module(id, surfaces);
-----------------
February 28, 2023
-----------------
PhreeqcRM: Revised names for PhreeqcRM test case source and output
files (Tests subdirectory of distribution). Added tests SimpleAdvect_cpp,
SimpleAdvect_c and SimpleAdvect_f90. All transport results are the same for
Advect_cpp, Advect_c, and Advect_f90, SimpleAdvect_cpp,
SimpleAdvect_c and SimpleAdvect_f90; however, the SimpleAdvect cases use
a minimal set of method calls, whereas the other cases demonstrate most
of the features of PhreeqcRM.
-----------------
February 26, 2023
-----------------
PhreeqcRM: Added method InitializeYAML to initialize a PhreeqcRM
instance.
It is possible to use a YAML file to initialize the PhreeqcRM instance.
If a GUI is used to define model input, a YAML file could be used to
transfer the data from the GUI to a PhreeqcRM instance. The YAML file
contains names of PhreeqcRM methods and associated data, for example:
RunFile:
workers: true
initial_phreeqc: true
utility: true
chemistry_name: advect.pqi
InitializeYAML can be used to process the directives defined in
the YAML file. The method InitializeYAML is equivalent to BMI_Initialize.
A C++ class named YAMLPhreeqcRM is included in the Tests subdirectory
of the PhreeqcRM distribution, which provides the capability to generate
a PhreeqcRM YAML file. The file WriteYAMLFile.cpp generates a YAML file
and advection_bmi_cpp.cpp reads and processes the file to initialize
a PhreeqcRM instance.
-----------------
February 26, 2023
-----------------
PhreeqcRM: Added a BMI (Basic Model Interface) for C++ and Fortran.
The interface is a repackaging of the available methods of
PhreeqcRM. All PhreeqcRM methods are available in addition
to the BMI methods.
New capabilities include (1) the capability to
retrieve units for the variables for BMI_GetValues and BMI_SetValues,
(2) the capability to use YAML (YAML ain't Markup Language)
to initialize a PhreeqcRM instance with the method BMI_Initialize
(which is equivalent to the method InitializeYAML), and (3) the availability
of pointers that always point to current variable values.
The YAML capability would be especially useful if a GUI (Graphical User Interface)
is used to set up model initial conditions. The GUI could write a YAML file
that contains directives for PhreeqcRM methods that need to be run and
the corresonding data needed to initialize a PhreeqcRM instance--for example,
setting units; running a PHREEQC input file to define initial and
boundary conditions; distribution of initial conditions to the model cells;
setting initial porosity, saturation, temperature, and pressure.
When a PhreeqcRM instance is created by the simulator, it can process
the YAML file with BMI_Initialize to execute the specified PhreeqcRM methods
to apply the data specified in the YAML file.
The following represents the way BMI methods would be used to implement
a sequential, noniterative transport calculation:
PhreeqcRM phreeqc_rm(nxyz, nthreads);
phreeqc_rm.Initialize("myfile.yaml");
int ncomps;
phreeqc_rm.GetValue("ComponentCount", &ncomps);
int ngrid;
phreeqc_rm.GetValue("GridCellCount", ngrid);
std::vector c(ngrid*ncomps, 0.0);
phreeqc_rm.GetValue("Concentrations", c.data());
phreeqc_rm.SetValue("TimeStep", 86400);
for(double time = 0; time < 864000; time+=86400)
{
// Take a transport time step here and update the vector c.
your_transport(c);
phreeqc_rm.SetValue("Time", time);
phreeqc_rm.SetValue("Concentrations", c.data());
phreeqc_rm.Update();
phreeqc_rm.GetValue("Concentrations", c.data());
}
The set of BMI methods is as follows:
std::string GetComponentName()
Returns "PhreeqcRM".
double GetCurrentTime()
Returns current time that has been set by the user.
double GetEndTime()
Returns current time plus the time step.
int GetInputItemCount()
Returns the number of variables that it is possible to set
with SetValue.
std::vector<std::string> GetInputVarNames()
Returns a list of the names of variables that can be set
with SetValue.
int GetOutputItemCount()
Returns the number of variables that it is possible to retrieve
with GetValue.
std::vector<std::string> GetOutputVarNames()
Returns a list of the names of variables that can be retrieved
with GetValue.
double GetTimeStep()
Returns the current time step that has been set by the user.
std::string GetTimeUnits()
Returns "seconds".
void GetValue(std::string name, void* dest)
Returns a value or vector of values for the variable identified by name
void GetValuePtr(std::string name, void* dest)
Returns a pointer to current values of a variable. This method
is available for selected variables.
int GetVarItemsize(std::string name)
Returns the number of bytes needed for one element of the variable
identified by name.
int GetVarNbytes(std::string name)
Returns the total number of bytes neded to store the value or vector
of values identified by name.
std::string GetVarType(std::string name)
Returns the type of the variable identified by name: "int", "double", or
"string".
std::string GetVarUnits(std::string name)
Returns the units associated with the variable identified by name.
void Initialize(std::string config_file)
Same as InitializeYAML discussed above.
void SetValue(std::string name, void* src)
Sets the value or vector of values for the variable identified by name.
void Update(void)
Calculates chemical reactions for a time step. It is equivalent to
the method RunCells. Equilibrium will be calculated between the solution
and all equilibrium reactants (EQUILIBRIUM_PHASES, EXCHANGE, etc), and
KINETICS will be integrated for the time step.
BMI is implemented in Fortran with USE BMIPhreeqcRM. Methods are nemed with a
prefix "bmif" and have an additional initial argument to identify the instance of
BMIPhreeqcRM that is being used.
-----------------
February 26, 2023
-----------------
PHREEQC: Fixed bug with Basic functions PR_P and PR_PHI. Values
were incorrect after the first step when INCREMENTAL_REACTIONS
was set to true.
-----------------
February 25, 2023
-----------------
ALL PROGRAMS: Added the latest version of the database Thermoddem to
the distributions of PHREEQC programs. The database was downloaded
from https://thermoddem.brgm.fr/.
-----------------
March 23, 2022
-----------------
PHREEQC: "MacInnes" was misspelled in one of the warning messages.
-----------------
December 18, 2021
-----------------
PHREEQC: Fixed transport bug where the end cell should not have
been processed in part of the calculation.
Full Changelog: v3.7.3...v3.8.0