diff --git a/phreeqc3-doc/.gitlab-ci.yml b/phreeqc3-doc/.gitlab-ci.yml index 45ffff7b..8f3a3913 100644 --- a/phreeqc3-doc/.gitlab-ci.yml +++ b/phreeqc3-doc/.gitlab-ci.yml @@ -41,12 +41,16 @@ trigger-downstream: script: - echo triggering iphreeqc - curl -X POST -F token=${IPHREEQC_TRIGGER} -F ref=master https://code.chs.usgs.gov/api/v4/projects/${IPHREEQC_ID}/trigger/pipeline + - sleep 180 - echo triggering iphreeqccom - curl -X POST -F token=${IPHREEQCCOM_TRIGGER} -F ref=master https://code.chs.usgs.gov/api/v4/projects/${IPHREEQCCOM_ID}/trigger/pipeline + - sleep 120 - echo triggering phast3-doc PHAST3_DOC_TRIGGER PHAST3_DOC_ID - curl -X POST -F token=${PHAST3_DOC_TRIGGER} -F ref=master https://code.chs.usgs.gov/api/v4/projects/${PHAST3_DOC_ID}/trigger/pipeline + - sleep 180 - echo triggering phreeqc3 PHREEQC3_TRIGGER PHREEQC3_ID - curl -X POST -F token=${PHREEQC3_TRIGGER} -F ref=master https://code.chs.usgs.gov/api/v4/projects/${PHREEQC3_ID}/trigger/pipeline + - sleep 360 ## Upstream Projects ## none diff --git a/phreeqc3-doc/RELEASE.TXT b/phreeqc3-doc/RELEASE.TXT index 391134bc..880a075b 100644 --- a/phreeqc3-doc/RELEASE.TXT +++ b/phreeqc3-doc/RELEASE.TXT @@ -1,4 +1,211 @@ Version @PHREEQC_VER@: @PHREEQC_DATE@ + + ----------------- + 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 smaller 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)22- species in + PHREEQC.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)22-). 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). + + Removed the NaCO3- species in PHREEQC.dat since they are not necessary for the calculation of + the specific conductance (SC) and their origin is unknown. Defined parameters in the + -analytical_expression, -gamma, -dw, -Vm and -viscosity for the NaHCO3 species in PHREEQC.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 CO32- 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 parameter. + 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 in c:\phreeqc\viscosity. + + 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: + Renamed GetDensity and related functions to GetDensityCalculated. + Renamed SetDensity and related functions to SetDensityUser. + + 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. + + Renamed GetSaturation and related functions to GetSaturationCalculated. + Renamed SetSaturation and related functions to SetSaturationUser. + + The values specified by SetSaturation 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& 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& 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& 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& 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); + ----------------- April 3, 2023 ----------------- @@ -59,9 +266,9 @@ Version @PHREEQC_VER@: @PHREEQC_DATE@ ----------------- February 26, 2023 ----------------- - PhreeqcRM: Added a BMI (Basic Model Interface) for C++. + 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 + PhreeqcRM. All PhreeqcRM methods are available in addition to the BMI methods. New capabilities include (1) the capability to @@ -82,95 +289,99 @@ Version @PHREEQC_VER@: @PHREEQC_DATE@ the YAML file with BMI_Initialize to execute the specified PhreeqcRM methods to apply the data specified in the YAML file. - The following is represents the way BMI methods would be used to implement + The following represents the way BMI methods would be used to implement a sequential, noniterative transport calculation: PhreeqcRM phreeqc_rm(nxyz, nthreads); - phreeqc_rm.BMI_Initialize("myfile.yaml"); + phreeqc_rm.Initialize("myfile.yaml"); int ncomps; - phreeqc_rm.BMI_GetValue("ComponentCount", &ncomps); + phreeqc_rm.GetValue("ComponentCount", &ncomps); int ngrid; - phreeqc_rm.BMI_GetValue("GridCellCount", ngrid); + phreeqc_rm.GetValue("GridCellCount", ngrid); std::vector c(ngrid*ncomps, 0.0); - phreeqc_rm.BMI_GetValue("Concentrations", c.data()); - phreeqc_rm.BMI_SetValue("TimeStep", 86400); + 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.BMI_SetValue("Time", time); - phreeqc_rm.BMI_SetValue("Concentrations", c.data()); - phreeqc_rm.BMI_Update(); - phreeqc_rm.BMI_GetValue("Concentrations", c.data()); + 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 BMI_GetComponentName() + std::string GetComponentName() Returns "PhreeqcRM". - double BMI_GetCurrentTime() + double GetCurrentTime() Returns current time that has been set by the user. - double BMI_GetEndTime() + double GetEndTime() Returns current time plus the time step. - int BMI_GetInputItemCount() + int GetInputItemCount() Returns the number of variables that it is possible to set - with BMI_SetValue. + with SetValue. - std::vector BMI_GetInputVarNames() + std::vector GetInputVarNames() Returns a list of the names of variables that can be set - with BMI_SetValue. + with SetValue. - int BMI_GetOutputItemCount() + int GetOutputItemCount() Returns the number of variables that it is possible to retrieve - with BMI_GetValue. + with GetValue. - std::vector BMI_GetOutputVarNames() + std::vector GetOutputVarNames() Returns a list of the names of variables that can be retrieved - with BMI_GetValue. + with GetValue. - double BMI_GetTimeStep() + double GetTimeStep() Returns the current time step that has been set by the user. - std::string BMI_GetTimeUnits() + std::string GetTimeUnits() Returns "seconds". - void BMI_GetValue(std::string name, void* dest) + void GetValue(std::string name, void* dest) Returns a value or vector of values for the variable identified by name - void BMI_GetValuePtr(std::string name, void* dest) + void GetValuePtr(std::string name, void* dest) Returns a pointer to current values of a variable. This method is available for selected variables. - int BMI_GetVarItemsize(std::string name) + int GetVarItemsize(std::string name) Returns the number of bytes needed for one element of the variable identified by name. - int BMI_GetVarNbytes(std::string 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 BMI_GetVarType(std::string name) + std::string GetVarType(std::string name) Returns the type of the variable identified by name: "int", "double", or "string". - std::string BMI_GetVarUnits(std::string name) + std::string GetVarUnits(std::string name) Returns the units associated with the variable identified by name. - void BMI_Initialize(std::string config_file) + void Initialize(std::string config_file) Same as InitializeYAML discussed above. - void BMI_SetValue(std::string name, void* src) + void SetValue(std::string name, void* src) Sets the value or vector of values for the variable identified by name. - void BMI_Update(void) + 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 -----------------