From 324810368e4c6393e3707f216112b9a8c2d8898f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Sep 2024 09:10:33 +0200 Subject: [PATCH 1/6] feat: remove header dependency to c++20 features --- poet/include/PhreeqcEngine.hpp | 76 ++----- poet/src/Engine.cpp | 200 ++++++++---------- poet/src/InitGetSet.cpp | 10 +- .../{ => Wrapper}/EquilibriumCompWrapper.cpp | 0 poet/src/{ => Wrapper}/EquilibriumWrapper.cpp | 0 .../Wrapper}/EquilibriumWrapper.hpp | 0 .../src/{ => Wrapper}/ExchangeCompWrapper.cpp | 0 poet/src/{ => Wrapper}/ExchangeWrapper.cpp | 0 .../Wrapper}/ExchangeWrapper.hpp | 0 poet/src/{ => Wrapper}/KineticCompWrapper.cpp | 0 poet/src/{ => Wrapper}/KineticWrapper.cpp | 0 .../Wrapper}/KineticWrapper.hpp | 0 poet/src/{ => Wrapper}/SolutionWrapper.cpp | 0 .../Wrapper}/SolutionWrapper.hpp | 0 .../{ => Wrapper}/SurfaceChargeWrapper.cpp | 0 poet/src/{ => Wrapper}/SurfaceCompWrapper.cpp | 0 poet/src/{ => Wrapper}/SurfaceWrapper.cpp | 0 .../Wrapper}/SurfaceWrapper.hpp | 0 poet/{include => src/Wrapper}/WrapperBase.hpp | 0 19 files changed, 103 insertions(+), 183 deletions(-) rename poet/src/{ => Wrapper}/EquilibriumCompWrapper.cpp (100%) rename poet/src/{ => Wrapper}/EquilibriumWrapper.cpp (100%) rename poet/{include => src/Wrapper}/EquilibriumWrapper.hpp (100%) rename poet/src/{ => Wrapper}/ExchangeCompWrapper.cpp (100%) rename poet/src/{ => Wrapper}/ExchangeWrapper.cpp (100%) rename poet/{include => src/Wrapper}/ExchangeWrapper.hpp (100%) rename poet/src/{ => Wrapper}/KineticCompWrapper.cpp (100%) rename poet/src/{ => Wrapper}/KineticWrapper.cpp (100%) rename poet/{include => src/Wrapper}/KineticWrapper.hpp (100%) rename poet/src/{ => Wrapper}/SolutionWrapper.cpp (100%) rename poet/{include => src/Wrapper}/SolutionWrapper.hpp (100%) rename poet/src/{ => Wrapper}/SurfaceChargeWrapper.cpp (100%) rename poet/src/{ => Wrapper}/SurfaceCompWrapper.cpp (100%) rename poet/src/{ => Wrapper}/SurfaceWrapper.cpp (100%) rename poet/{include => src/Wrapper}/SurfaceWrapper.hpp (100%) rename poet/{include => src/Wrapper}/WrapperBase.hpp (100%) diff --git a/poet/include/PhreeqcEngine.hpp b/poet/include/PhreeqcEngine.hpp index fd455d90..91ae9a2c 100644 --- a/poet/include/PhreeqcEngine.hpp +++ b/poet/include/PhreeqcEngine.hpp @@ -1,49 +1,35 @@ - #pragma once -#include "EquilibriumWrapper.hpp" -#include "ExchangeWrapper.hpp" -#include "KineticWrapper.hpp" #include "POETInit.hpp" -#include "SolutionWrapper.hpp" -#include "SurfaceWrapper.hpp" -#include - -#include -#include -#include -#include -#include -#include -#include #include -#include -#include #include /** * @brief Class for running Phreeqc wrappped in POET * - * Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts to + * Direct interface to Phreeqc, without utilizing *_MODIFY keywords/scripts + to * set new values. Use with already initialized Phreeqc config. * */ -class PhreeqcEngine : public IPhreeqc { +class PhreeqcEngine { public: /** * @brief Construct a new Phreeqc Engine object * - * Construct a new Phreeqc Engine object by previously initialized POETConfig. + * Construct a new Phreeqc Engine object by previously initialized + POETConfig. * * @param config Holds the configuration for the Phreeqc engine. */ - PhreeqcEngine(const POETConfig &config) : IPhreeqc() { - this->LoadDatabaseString(config.database.c_str()); - this->RunString(config.input_script.c_str()); + PhreeqcEngine(const POETConfig &config); - this->init_wrappers(config.cell); - } + /** + * @brief Destroy the Phreeqc Engine object + * + */ + ~PhreeqcEngine(); /** * @brief Siimulate a cell for a given time step @@ -56,42 +42,6 @@ public: void runCell(std::vector &cell_values, double time_step); private: - void init_wrappers(const POETInitCell &cell); - void run(double time_step); - - cxxSolution *Get_solution(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n); - } - - cxxExchange *Get_exchange(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n); - } - - cxxKinetics *Get_kinetic(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n); - } - - cxxPPassemblage *Get_equilibrium(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(), - n); - } - - cxxSurface *Get_surface(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n); - } - - void get_essential_values(std::span &data); - - void set_essential_values(const std::span &data); - - std::unique_ptr solutionWrapperPtr; - std::unique_ptr exchangeWrapperPtr; - std::unique_ptr kineticsWrapperPtr; - std::unique_ptr equilibriumWrapperPtr; - std::unique_ptr surfaceWrapperPtr; - - bool has_exchange = false; - bool has_kinetics = false; - bool has_equilibrium = false; - bool has_surface = false; + class Impl; + std::unique_ptr impl; }; \ No newline at end of file diff --git a/poet/src/Engine.cpp b/poet/src/Engine.cpp index f3af6f1d..c74b73da 100644 --- a/poet/src/Engine.cpp +++ b/poet/src/Engine.cpp @@ -4,7 +4,89 @@ #include #include -void PhreeqcEngine::init_wrappers(const POETInitCell &cell) { +#include +#include + +#include "Wrapper/EquilibriumWrapper.hpp" +#include "Wrapper/ExchangeWrapper.hpp" +#include "Wrapper/KineticWrapper.hpp" +#include "Wrapper/SolutionWrapper.hpp" +#include "Wrapper/SurfaceWrapper.hpp" +class PhreeqcEngine::Impl : public IPhreeqc { +public: + Impl(const POETConfig &config); + void init_wrappers(const POETInitCell &cell); + void run(double time_step); + + cxxSolution *Get_solution(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n); + } + + cxxExchange *Get_exchange(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n); + } + + cxxKinetics *Get_kinetic(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n); + } + + cxxPPassemblage *Get_equilibrium(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(), + n); + } + + cxxSurface *Get_surface(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n); + } + + void get_essential_values(std::span &data); + + void set_essential_values(const std::span &data); + + std::unique_ptr solutionWrapperPtr; + std::unique_ptr exchangeWrapperPtr; + std::unique_ptr kineticsWrapperPtr; + std::unique_ptr equilibriumWrapperPtr; + std::unique_ptr surfaceWrapperPtr; + + bool has_exchange = false; + bool has_kinetics = false; + bool has_equilibrium = false; + bool has_surface = false; +}; + +PhreeqcEngine::PhreeqcEngine(const POETConfig &config) + : impl(std::make_unique(config)) {} + +PhreeqcEngine::Impl::Impl(const POETConfig &config) { + this->LoadDatabaseString(config.database.c_str()); + this->RunString(config.input_script.c_str()); + + this->init_wrappers(config.cell); +} + +PhreeqcEngine::~PhreeqcEngine() = default; + +void PhreeqcEngine::runCell(std::vector &cell_values, + double time_step) { + // skip ID + std::span cell_data{cell_values.begin() + 1, cell_values.end()}; + + this->impl->set_essential_values(cell_data); + this->impl->run(time_step); + this->impl->get_essential_values(cell_data); +} + +void PhreeqcEngine::Impl::run(double time_step) { + std::stringstream time_ss; + time_ss << std::fixed << std::setprecision(20) << time_step; + + const std::string runs_string = + "RUN_CELLS\n -cells 1\n -time_step " + time_ss.str() + "\nEND\n"; + this->RunString(runs_string.c_str()); +} + +void PhreeqcEngine::Impl::init_wrappers(const POETInitCell &cell) { // Solutions this->solutionWrapperPtr = @@ -41,7 +123,7 @@ void PhreeqcEngine::init_wrappers(const POETInitCell &cell) { } } -void PhreeqcEngine::get_essential_values(std::span &data) { +void PhreeqcEngine::Impl::get_essential_values(std::span &data) { this->solutionWrapperPtr->get(data); @@ -76,54 +158,9 @@ void PhreeqcEngine::get_essential_values(std::span &data) { data.subspan(offset, this->surfaceWrapperPtr->size())}; this->surfaceWrapperPtr->get(surf_span); } - - // // Exchange - // if (this->Get_exchange(cell_number) != NULL) { - // std::vector exch_values(this->exchangeWrapperPtr->size()); - - // std::span exch_span{exch_values}; - - // this->exchangeWrapperPtr->get(exch_span); - - // // this->Get_exchange(cell_number)->get_essential_values(exch_values); - // essentials.insert(essentials.end(), exch_values.begin(), - // exch_values.end()); - // } - - // // Kinetics - // if (this->Get_kinetic(cell_number) != NULL) { - // std::vector kin_values; - - // this->Get_kinetic(cell_number)->get_essential_values(kin_values); - // essentials.insert(essentials.end(), kin_values.begin(), - // kin_values.end()); - // } - - // // PPassemblage - // if (this->Get_equilibrium(cell_number) != NULL) { - // std::vector equ_values; - - // this->Get_equilibrium(cell_number)->get_essential_values(equ_values); - // essentials.insert(essentials.end(), equ_values.begin(), - // equ_values.end()); - // } - - // // Surface - // if (this->Get_surface(cell_number) != NULL) { - // std::vector surf_values(this->surfaceWrapperPtr->size()); - // std::span surf_span{surf_values}; - - // this->surfaceWrapperPtr->get(surf_span); - - // // this->Get_surface(cell_number)->get_essential_values(surf_values); - // essentials.insert(essentials.end(), surf_values.begin(), - // surf_values.end()); - // } - - // return essentials; } -void PhreeqcEngine::set_essential_values(const std::span &data) { +void PhreeqcEngine::Impl::set_essential_values(const std::span &data) { this->solutionWrapperPtr->set(data); this->PhreeqcPtr->initial_solutions_poet(1); @@ -159,71 +196,4 @@ void PhreeqcEngine::set_essential_values(const std::span &data) { data.subspan(offset, this->surfaceWrapperPtr->size())}; this->surfaceWrapperPtr->set(surf_span); } - - // // Solutions - // if (this->Get_solution(cell_number) != NULL) { - // this->Get_solution(cell_number)->set_essential_values(dump_it, order); - // this->PhreeqcPtr->initial_solutions_poet(cell_number); - // } - - // // Exchange - // if (this->Get_exchange(cell_number) != NULL) { - // auto &exch_pointer = this->exchangeWrapperPtr; - // std::span exch_span{dump_it, exch_pointer->size()}; - // exch_pointer->set(exch_span); - - // dump_it += exch_pointer->size(); - - // // this->Get_exchange(cell_number)->set_essential_values(dump_it); - // // this->PhreeqcPtr->Rxn_new_exchange.insert(cell_number); - // // this->Get_exchange(cell_number)->Set_new_def(true); - // // this->Get_exchange(cell_number)->Set_solution_equilibria(true); - // // this->Get_exchange(cell_number)->Set_n_solution(cell_number); - // // this->PhreeqcPtr->initial_exchangers(FALSE); - // } - - // // Kinetics - // if (this->Get_kinetic(cell_number) != NULL) { - // this->Get_kinetic(cell_number)->set_essential_values(dump_it); - // } - - // // PPassemblage - // if (this->Get_equilibrium(cell_number) != NULL) { - // this->Get_equilibrium(cell_number)->set_essential_values(dump_it); - // } - - // // Surface - // if (this->Get_surface(cell_number) != NULL) { - // auto *tmp = this->Get_surface(cell_number); - // auto &surf_pointer = this->surfaceWrapperPtr; - // std::span surf_span{dump_it, surf_pointer->size()}; - // surf_pointer->set(surf_span); - - // dump_it += surf_pointer->size(); - - // // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number); - - // // this->Get_surface(cell_number)->set_essential_values(dump_it); - // // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number); - // // this->PhreeqcPtr->initial_surfaces(TRUE); - // } -} - -void PhreeqcEngine::runCell(std::vector &cell_values, - double time_step) { - // skip ID - std::span cell_data{cell_values.begin() + 1, cell_values.end()}; - - this->set_essential_values(cell_data); - this->run(time_step); - this->get_essential_values(cell_data); -} - -void PhreeqcEngine::run(double time_step) { - std::stringstream time_ss; - time_ss << std::fixed << std::setprecision(20) << time_step; - - const std::string runs_string = - "RUN_CELLS\n -cells 1\n -time_step " + time_ss.str() + "\nEND\n"; - this->RunString(runs_string.c_str()); } \ No newline at end of file diff --git a/poet/src/InitGetSet.cpp b/poet/src/InitGetSet.cpp index 04bb03f3..1f5a813a 100644 --- a/poet/src/InitGetSet.cpp +++ b/poet/src/InitGetSet.cpp @@ -1,9 +1,9 @@ -#include "EquilibriumWrapper.hpp" -#include "ExchangeWrapper.hpp" -#include "KineticWrapper.hpp" #include "PhreeqcInit.hpp" -#include "SolutionWrapper.hpp" -#include "SurfaceWrapper.hpp" +#include "Wrapper/EquilibriumWrapper.hpp" +#include "Wrapper/ExchangeWrapper.hpp" +#include "Wrapper/KineticWrapper.hpp" +#include "Wrapper/SolutionWrapper.hpp" +#include "Wrapper/SurfaceWrapper.hpp" #include "global_structures.h" #include diff --git a/poet/src/EquilibriumCompWrapper.cpp b/poet/src/Wrapper/EquilibriumCompWrapper.cpp similarity index 100% rename from poet/src/EquilibriumCompWrapper.cpp rename to poet/src/Wrapper/EquilibriumCompWrapper.cpp diff --git a/poet/src/EquilibriumWrapper.cpp b/poet/src/Wrapper/EquilibriumWrapper.cpp similarity index 100% rename from poet/src/EquilibriumWrapper.cpp rename to poet/src/Wrapper/EquilibriumWrapper.cpp diff --git a/poet/include/EquilibriumWrapper.hpp b/poet/src/Wrapper/EquilibriumWrapper.hpp similarity index 100% rename from poet/include/EquilibriumWrapper.hpp rename to poet/src/Wrapper/EquilibriumWrapper.hpp diff --git a/poet/src/ExchangeCompWrapper.cpp b/poet/src/Wrapper/ExchangeCompWrapper.cpp similarity index 100% rename from poet/src/ExchangeCompWrapper.cpp rename to poet/src/Wrapper/ExchangeCompWrapper.cpp diff --git a/poet/src/ExchangeWrapper.cpp b/poet/src/Wrapper/ExchangeWrapper.cpp similarity index 100% rename from poet/src/ExchangeWrapper.cpp rename to poet/src/Wrapper/ExchangeWrapper.cpp diff --git a/poet/include/ExchangeWrapper.hpp b/poet/src/Wrapper/ExchangeWrapper.hpp similarity index 100% rename from poet/include/ExchangeWrapper.hpp rename to poet/src/Wrapper/ExchangeWrapper.hpp diff --git a/poet/src/KineticCompWrapper.cpp b/poet/src/Wrapper/KineticCompWrapper.cpp similarity index 100% rename from poet/src/KineticCompWrapper.cpp rename to poet/src/Wrapper/KineticCompWrapper.cpp diff --git a/poet/src/KineticWrapper.cpp b/poet/src/Wrapper/KineticWrapper.cpp similarity index 100% rename from poet/src/KineticWrapper.cpp rename to poet/src/Wrapper/KineticWrapper.cpp diff --git a/poet/include/KineticWrapper.hpp b/poet/src/Wrapper/KineticWrapper.hpp similarity index 100% rename from poet/include/KineticWrapper.hpp rename to poet/src/Wrapper/KineticWrapper.hpp diff --git a/poet/src/SolutionWrapper.cpp b/poet/src/Wrapper/SolutionWrapper.cpp similarity index 100% rename from poet/src/SolutionWrapper.cpp rename to poet/src/Wrapper/SolutionWrapper.cpp diff --git a/poet/include/SolutionWrapper.hpp b/poet/src/Wrapper/SolutionWrapper.hpp similarity index 100% rename from poet/include/SolutionWrapper.hpp rename to poet/src/Wrapper/SolutionWrapper.hpp diff --git a/poet/src/SurfaceChargeWrapper.cpp b/poet/src/Wrapper/SurfaceChargeWrapper.cpp similarity index 100% rename from poet/src/SurfaceChargeWrapper.cpp rename to poet/src/Wrapper/SurfaceChargeWrapper.cpp diff --git a/poet/src/SurfaceCompWrapper.cpp b/poet/src/Wrapper/SurfaceCompWrapper.cpp similarity index 100% rename from poet/src/SurfaceCompWrapper.cpp rename to poet/src/Wrapper/SurfaceCompWrapper.cpp diff --git a/poet/src/SurfaceWrapper.cpp b/poet/src/Wrapper/SurfaceWrapper.cpp similarity index 100% rename from poet/src/SurfaceWrapper.cpp rename to poet/src/Wrapper/SurfaceWrapper.cpp diff --git a/poet/include/SurfaceWrapper.hpp b/poet/src/Wrapper/SurfaceWrapper.hpp similarity index 100% rename from poet/include/SurfaceWrapper.hpp rename to poet/src/Wrapper/SurfaceWrapper.hpp diff --git a/poet/include/WrapperBase.hpp b/poet/src/Wrapper/WrapperBase.hpp similarity index 100% rename from poet/include/WrapperBase.hpp rename to poet/src/Wrapper/WrapperBase.hpp From 6fd2e96411853bf7d43b343c8c167308d28d6ddd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Sep 2024 09:11:28 +0200 Subject: [PATCH 2/6] refactor: Update C++ standard to c++20 in IPhreeqcPOET target --- poet/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/poet/CMakeLists.txt b/poet/CMakeLists.txt index f2125ba1..f6f4b0d6 100644 --- a/poet/CMakeLists.txt +++ b/poet/CMakeLists.txt @@ -9,4 +9,4 @@ target_include_directories(IPhreeqcPOET PUBLIC PRIVATE src) # Set C++23 standard -target_compile_features(IPhreeqcPOET PUBLIC cxx_std_23) \ No newline at end of file +target_compile_features(IPhreeqcPOET PUBLIC cxx_std_20) \ No newline at end of file From e42c120518ad1f69a421e8b18352bd97ae168651 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Sep 2024 15:08:55 +0200 Subject: [PATCH 3/6] Add 'app' directory from `golem` branch --- CMakeLists.txt | 1 + app/CMakeLists.txt | 25 + app/ipqEngine.cpp | 122 +++ app/ipqMain.cpp | 92 +++ app/phreeqc.dat | 1853 ++++++++++++++++++++++++++++++++++++++++++++ app/testEngine.pqi | 13 + 6 files changed, 2106 insertions(+) create mode 100644 app/CMakeLists.txt create mode 100644 app/ipqEngine.cpp create mode 100644 app/ipqMain.cpp create mode 100644 app/phreeqc.dat create mode 100644 app/testEngine.pqi diff --git a/CMakeLists.txt b/CMakeLists.txt index 0fa431a5..5f4bb41a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -393,6 +393,7 @@ if (BUILD_CLR_LIBS) endif() add_subdirectory(poet) +add_subdirectory(app) include (CTest) diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt new file mode 100644 index 00000000..8055e77d --- /dev/null +++ b/app/CMakeLists.txt @@ -0,0 +1,25 @@ +## Time-stamp: "Last modified 2024-08-26 17:42:18 delucia" + +add_executable(ipqMain + ipqMain.cpp + ) +add_executable(ipqEngine + ipqEngine.cpp + ) + +# add_library(IPhreeqcPOET ${PQC_POET_SRC}) +##target_link_libraries(IPhreeqcPOET PUBLIC IPhreeqc) + +target_include_directories(IPhreeqcPOET PUBLIC + $ + $ + PRIVATE src) + + +target_link_libraries(ipqMain + IPhreeqcPOET +) + +target_link_libraries(ipqEngine + IPhreeqcPOET +) diff --git a/app/ipqEngine.cpp b/app/ipqEngine.cpp new file mode 100644 index 00000000..7199ea92 --- /dev/null +++ b/app/ipqEngine.cpp @@ -0,0 +1,122 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "POETInit.hpp" +#include "PhreeqcEngine.hpp" +#include "PhreeqcInit.hpp" + +static std::string readFile(const std::string &path) { + std::string string_rpath(PATH_MAX, '\0'); + + if (realpath(path.c_str(), string_rpath.data()) == nullptr) { + throw std::runtime_error("Failed to resolve the realpath to file " + path); + } + + std::ifstream file(string_rpath); + + if (!file.is_open()) { + throw std::runtime_error("Failed to open file: " + path); + } + + std::stringstream buffer; + buffer << file.rdbuf(); + + return buffer.str(); +} + +void printVectorsAsTable(const std::vector &names, + const std::vector &values) { + if (names.size() != values.size()) { + std::cerr << "Error: Vectors are not of the same length." << std::endl; + return; + } + + std::cout << std::left << std::setw(20) << "Name" << std::setw(20) << "Value" + << std::endl; + std::cout << std::string(40, '-') << std::endl; + + for (std::size_t i = 0; i < names.size(); ++i) { + std::cout << std::left << std::setw(20) << names[i] << std::setw(20) + << values[i] << std::endl; + } +} + +template void printVector(const std::vector &vec) { + for (const auto &elem : vec) { + std::cout << elem << " "; + } + std::cout << std::endl; +} + +int main(int argc, char *argv[]) { + + if (argc < 3) { + std::cout << "Two args needed, PHREEQC script and database" << std::endl; + return 1; + } + + const std::string &script = argv[1]; + const std::string &database = argv[2]; + wordexp_t exp_result; + + std::cout << "Reading script from file " << script << " using database " + << database << std::endl; + + // expand the script and database path + const std::string read_script = readFile(script); + const std::string read_database = readFile(database); + + // initialize the module instance "ipqcmod" + PhreeqcInit ipqcmod(read_database, read_script); + + // allocate tabular data + PhreeqcInit::PhreeqcMat pqc_mat = ipqcmod.getPhreeqcMat(); + + // for (std::size_t i = 0; i < nrows; ++i) { + auto pure_names = ipqcmod.getEquilibriumNames(pqc_mat.ids[0]); + auto totals_names = ipqcmod.getSolutionNames(pqc_mat.ids[0]); + + POETConfig MyConfig; + MyConfig.database = read_database; + MyConfig.input_script = read_script; + MyConfig.cell = POETInitCell{ + // H, O and charge need to be removed. I'm not sure why I did it this way. + std::vector(totals_names.begin() + 3, totals_names.end()), + {}, + {}, + {}, + pure_names, + {}, + {}}; + + std::vector full_names = {"ID"}; + full_names.insert(full_names.end(), pqc_mat.names.begin(), + pqc_mat.names.end()); + + PhreeqcEngine MyEngine(MyConfig); + std::vector NewInp = pqc_mat.values[0]; + NewInp.insert(NewInp.begin(), 1); + + std::cout << "Before running the cell:\n"; + printVectorsAsTable(full_names, NewInp); + + // std::for_each(NewInp.begin()+3, NewInp.end(), + // [](double a) { return a * 0.8;}); + + NewInp[7] += 0.2; // Cl + NewInp[9] += 0.2; // Na + MyEngine.runCell(NewInp, 200); + + std::cout << "\n\nAfter running the cell:\n"; + printVectorsAsTable(full_names, NewInp); + + return 0; +} diff --git a/app/ipqMain.cpp b/app/ipqMain.cpp new file mode 100644 index 00000000..f38b68a8 --- /dev/null +++ b/app/ipqMain.cpp @@ -0,0 +1,92 @@ +// Time-stamp: "Last modified 2024-08-27 09:00:18 delucia" +#include +#include +#include +#include +#include +#include +#include + +#include "POETInit.hpp" +#include "PhreeqcEngine.hpp" +#include "PhreeqcInit.hpp" + +static std::string readFile(const std::string &path) { + std::string string_rpath(PATH_MAX, '\0'); + + if (realpath(path.c_str(), string_rpath.data()) == nullptr) { + throw std::runtime_error("Failed to resolve the realpath to file " + path); + } + + std::ifstream file(string_rpath); + + if (!file.is_open()) { + throw std::runtime_error("Failed to open file: " + path); + } + + std::stringstream buffer; + buffer << file.rdbuf(); + + return buffer.str(); +} + +int main(int argc, char *argv[]) { + + const std::string &script = argv[1]; + const std::string &database = argv[2]; + // wordexp_t exp_result; + + std::cout << "Reading script from file " << script << " using database " + << database << std::endl; + + // expand the script and database path + const std::string read_script = readFile(script); + const std::string read_database = readFile(database); + + // initialize the module instance "ipqcmod" + PhreeqcInit ipqcmod(read_database, read_script); + + // allocate tabular data + PhreeqcInit::PhreeqcMat pqc_mat = ipqcmod.getPhreeqcMat(); + + const std::size_t ncols = pqc_mat.names.size(); + const std::size_t nrows = pqc_mat.values.size(); + + std::cout << "Script contains " << nrows << " rows and " << ncols + << " columns" << std::endl; + + // pqc_mat.names.insert(pqc_mat.names.begin(), "ID"); + + // Rcpp::NumericMatrix mat(nrows, ncols + 1); + + std::cout << std::endl << "Number of CELLS/SOLUTIONS: " << nrows << std::endl; + + for (std::size_t i = 0; i < nrows; ++i) { + auto pure_names = ipqcmod.getEquilibriumNames(pqc_mat.ids[i]); + auto totals_names = ipqcmod.getSolutionNames(pqc_mat.ids[i]); + + std::cout << std::endl << "SOL " << pqc_mat.ids[i] << " - Pure phases : "; + for (std::size_t j = 0; j < pure_names.size(); j++) { + std::cout << pure_names[j] << ", "; + } + + std::cout << std::endl << "SOL " << pqc_mat.ids[i] << " - Total concs: "; + for (std::size_t j = 0; j < totals_names.size(); j++) { + std::cout << totals_names[j] << ", "; + } + std::cout << std::endl; + } + + // alternatively, use the pqc_mat object + std::cout << std::endl << "Total pqc_mat: " << std::endl; + for (std::size_t i = 0; i < nrows; ++i) { + std::cout << "ID: " << pqc_mat.ids[i] << " - "; + for (std::size_t j = 0; j < ncols; ++j) { + std::cout << pqc_mat.names[j] << ": " << pqc_mat.values[i][j] << ", "; + } + std::cout << std::endl; + } + std::cout << std::endl; + + return 0; +} diff --git a/app/phreeqc.dat b/app/phreeqc.dat new file mode 100644 index 00000000..9e69a187 --- /dev/null +++ b/app/phreeqc.dat @@ -0,0 +1,1853 @@ +# PHREEQC.DAT for calculating pressure dependence of reactions, with +# molal volumina of aqueous species and of minerals, and +# critical temperatures and pressures of gases used in Peng-Robinson's EOS. +# Details are given at the end of this file. + +SOLUTION_MASTER_SPECIES +# +#element species alk gfw_formula element_gfw +# +H H+ -1.0 H 1.008 +H(0) H2 0 H +H(1) H+ -1.0 0 +E e- 0 0.0 0 +O H2O 0 O 16.0 +O(0) O2 0 O +O(-2) H2O 0 0 +Ca Ca+2 0 Ca 40.08 +Mg Mg+2 0 Mg 24.312 +Na Na+ 0 Na 22.9898 +K K+ 0 K 39.102 +Fe Fe+2 0 Fe 55.847 +Fe(+2) Fe+2 0 Fe +Fe(+3) Fe+3 -2.0 Fe +Mn Mn+2 0 Mn 54.938 +Mn(+2) Mn+2 0 Mn +Mn(+3) Mn+3 0 Mn +Al Al+3 0 Al 26.9815 +Ba Ba+2 0 Ba 137.34 +Sr Sr+2 0 Sr 87.62 +Si H4SiO4 0 SiO2 28.0843 +Cl Cl- 0 Cl 35.453 +C CO3-2 2.0 HCO3 12.0111 +C(+4) CO3-2 2.0 HCO3 +C(-4) CH4 0 CH4 +Alkalinity CO3-2 1.0 Ca0.5(CO3)0.5 50.05 +S SO4-2 0 SO4 32.064 +S(6) SO4-2 0 SO4 +S(-2) HS- 1.0 S +N NO3- 0 N 14.0067 +N(+5) NO3- 0 N +N(+3) NO2- 0 N +N(0) N2 0 N +N(-3) NH4+ 0 N 14.0067 +#Amm AmmH+ 0 AmmH 17.031 +B H3BO3 0 B 10.81 +P PO4-3 2.0 P 30.9738 +F F- 0 F 18.9984 +Li Li+ 0 Li 6.939 +Br Br- 0 Br 79.904 +Zn Zn+2 0 Zn 65.37 +Cd Cd+2 0 Cd 112.4 +Pb Pb+2 0 Pb 207.19 +Cu Cu+2 0 Cu 63.546 +Cu(+2) Cu+2 0 Cu +Cu(+1) Cu+1 0 Cu +# redox-uncoupled gases +Hdg Hdg 0 Hdg 2.016 # H2 gas +Oxg Oxg 0 Oxg 32 # O2 gas +Mtg Mtg 0 Mtg 16.032 # CH4 gas +Sg H2Sg 1.0 H2Sg 34.08 +Ntg Ntg 0 Ntg 28.0134 # N2 gas + +SOLUTION_SPECIES +H+ = H+ + -gamma 9.0 0 + -dw 9.31e-9 1000 0.46 1e-10 # The dw parameters are defined in ref. 3. +# Dw(TK) = 9.31e-9 * exp(1000 / TK - 1000 / 298.15) * TK * 0.89 / (298.15 * viscos) +# Dw(I) = Dw(TK) * exp(-0.46 * DH_A * |z_H+| * I^0.5 / (1 + DH_B * I^0.5 * 1e-10 / (1 + I^0.75))) +e- = e- +H2O = H2O +# H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence +Ca+2 = Ca+2 + -gamma 5.0 0.1650 + -dw 0.793e-9 97 3.4 24.6 + -Vm -0.3456 -7.252 6.149 -2.479 1.239 5 1.60 -57.1 -6.12e-3 1 # ref. 1 +Mg+2 = Mg+2 + -gamma 5.5 0.20 + -dw 0.705e-9 111 2.4 13.7 + -Vm -1.410 -8.6 11.13 -2.39 1.332 5.5 1.29 -32.9 -5.86e-3 1 # ref. 1 +Na+ = Na+ + -gamma 4.0 0.075 + -gamma 4.08 0.082 # halite solubility + -dw 1.33e-9 122 1.52 3.70 + -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.566 # ref. 1 +# for calculating densities (rho) when I > 3... + # -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.45 +K+ = K+ + -gamma 3.5 0.015 + -dw 1.96e-9 395 2.5 21 + -Vm 3.322 -1.473 6.534 -2.712 9.06e-2 3.5 0 29.7 0 1 # ref. 1 +Fe+2 = Fe+2 + -gamma 6.0 0 + -dw 0.719e-9 + -Vm -0.3255 -9.687 1.536 -2.379 0.3033 6 -4.21e-2 39.7 0 1 # ref. 1 +Mn+2 = Mn+2 + -gamma 6.0 0 + -dw 0.688e-9 + -Vm -1.10 -8.03 4.08 -2.45 1.4 6 8.07 0 -1.51e-2 0.118 # ref. 2 +Al+3 = Al+3 + -gamma 9.0 0 + -dw 0.559e-9 + -Vm -2.28 -17.1 10.9 -2.07 2.87 9 0 0 5.5e-3 1 # ref. 2 and Barta and Hepler, 1986, Can. J.C. 64, 353. +Ba+2 = Ba+2 + -gamma 5.0 0 + -gamma 4.0 0.153 # Barite solubility + -dw 0.848e-9 46 + -Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -8.35e-3 1 # ref. 1 +Sr+2 = Sr+2 + -gamma 5.260 0.121 + -dw 0.794e-9 161 + -Vm -1.57e-2 -10.15 10.18 -2.36 0.860 5.26 0.859 -27.0 -4.1e-3 1.97 # ref. 1 +H4SiO4 = H4SiO4 + -dw 1.10e-9 + -Vm 10.5 1.7 20 -2.7 0.1291 # supcrt + 2*H2O in a1 +Cl- = Cl- + -gamma 3.5 0.015 + -gamma 3.63 0.017 # cf. pitzer.dat + -dw 2.03e-9 194 1.6 6.9 + -Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1 # ref. 1 +CO3-2 = CO3-2 + -gamma 5.4 0 + -dw 0.955e-9 0 1.12 2.84 + -Vm 5.95 0 0 -5.67 6.85 0 1.37 106 -0.0343 1 # ref. 1 +SO4-2 = SO4-2 + -gamma 5.0 -0.04 + -dw 1.07e-9 34 2.08 13.4 + -Vm 8.0 2.3 -46.04 6.245 3.82 0 0 0 0 1 # ref. 1 +NO3- = NO3- + -gamma 3.0 0 + -dw 1.9e-9 184 1.85 3.85 + -Vm 6.32 6.78 0 -3.06 0.346 0 0.93 0 -0.012 1 # ref. 1 +#AmmH+ = AmmH+ +# -gamma 2.5 0 +# -dw 1.98e-9 312 0.95 4.53 +# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # ref. 1 +H3BO3 = H3BO3 + -dw 1.1e-9 + -Vm 7.0643 8.8547 3.5844 -3.1451 -.2000 # supcrt +PO4-3 = PO4-3 + -gamma 4.0 0 + -dw 0.612e-9 + -Vm 1.24 -9.07 9.31 -2.4 5.61 0 0 0 -1.41e-2 1 # ref. 2 +F- = F- + -gamma 3.5 0 + -dw 1.46e-9 + -Vm 0.928 1.36 6.27 -2.84 1.84 0 0 -0.318 0 1 # ref. 2 +Li+ = Li+ + -gamma 6.0 0 + -dw 1.03e-9 80 + -Vm -0.419 -0.069 13.16 -2.78 0.416 0 0.296 -12.4 -2.74e-3 1.26 # ref. 2 and Ellis, 1968, J. Chem. Soc. A, 1138 +Br- = Br- + -gamma 3.0 0 + -dw 2.01e-9 258 + -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 # ref. 2 +Zn+2 = Zn+2 + -gamma 5.0 0 + -dw 0.715e-9 + -Vm -1.96 -10.4 14.3 -2.35 1.46 5 -1.43 24 1.67e-2 1.11 # ref. 2 +Cd+2 = Cd+2 + -dw 0.717e-9 + -Vm 1.63 -10.7 1.01 -2.34 1.47 5 0 0 0 1 # ref. 2 +Pb+2 = Pb+2 + -dw 0.945e-9 + -Vm -.0051 -7.7939 8.8134 -2.4568 1.0788 4.5 # supcrt +Cu+2 = Cu+2 + -gamma 6.0 0 + -dw 0.733e-9 + -Vm -1.13 -10.5 7.29 -2.35 1.61 6 9.78e-2 0 3.42e-3 1 # ref. 2 +# redox-uncoupled gases +Hdg = Hdg # H2 + -dw 5.13e-9 + -Vm 6.52 0.78 0.12 # supcrt +Oxg = Oxg # O2 + -dw 2.35e-9 + -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt +Mtg = Mtg # CH4 + -dw 1.85e-9 + -Vm 9.01 -1.11 0 -1.85 -1.50 # ref. 1 + Hnedkovsky et al., 1996, JCT 28, 125 +Ntg = Ntg # N2 + -dw 1.96e-9 + -Vm 7 # Pray et al., 1952, IEC 44. 1146 +H2Sg = H2Sg # H2S + -dw 2.1e-9 + -Vm 1.39 28.3 0 -7.22 -0.59 # ref. 1 + Hnedkovsky et al., 1996, JCT 28, 125 +# aqueous species +H2O = OH- + H+ + -analytic 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-5 + -gamma 3.5 0 + -dw 5.27e-9 548 0.52 1e-10 + -Vm -9.66 28.5 80.0 -22.9 1.89 0 1.09 0 0 1 # ref. 1 +2 H2O = O2 + 4 H+ + 4 e- + -log_k -86.08 + -delta_h 134.79 kcal + -dw 2.35e-9 + -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt +2 H+ + 2 e- = H2 + -log_k -3.15 + -delta_h -1.759 kcal + -dw 5.13e-9 + -Vm 6.52 0.78 0.12 # supcrt +CO3-2 + H+ = HCO3- + -log_k 10.329 + -delta_h -3.561 kcal + -analytic 107.8871 0.03252849 -5151.79 -38.92561 563713.9 + -gamma 5.4 0 + -dw 1.18e-9 0 1.43 1e-10 + -Vm 8.472 0 -11.5 0 1.56 0 0 146 3.16e-3 1 # ref. 1 +CO3-2 + 2 H+ = CO2 + H2O + -log_k 16.681 + -delta_h -5.738 kcal + -analytic 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 + -dw 1.92e-9 + -Vm 7.29 0.92 2.07 -1.23 -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171 +2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T + -log_k -1.8 + -analytical_expression 8.68 -0.0103 -2190 + -Vm 14.58 1.84 4.14 -2.46 -3.20 +CO3-2 + 10 H+ + 8 e- = CH4 + 3 H2O + -log_k 41.071 + -delta_h -61.039 kcal + -dw 1.85e-9 + -Vm 9.01 -1.11 0 -1.85 -1.50 # ref. 1 + Hnedkovsky et al., 1996, JCT 28, 125 +SO4-2 + H+ = HSO4- + -log_k 1.988 + -delta_h 3.85 kcal + -analytic -56.889 0.006473 2307.9 19.8858 + -dw 1.33e-9 + -Vm 8.2 9.2590 2.1108 -3.1618 1.1748 0 -0.3 15 0 1 # ref. 1 +HS- = S-2 + H+ + -log_k -12.918 + -delta_h 12.1 kcal + -gamma 5.0 0 + -dw 0.731e-9 +SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O + -log_k 33.65 + -delta_h -60.140 kcal + -gamma 3.5 0 + -dw 1.73e-9 + -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt +HS- + H+ = H2S + -log_k 6.994 + -delta_h -5.30 kcal + -analytical -11.17 0.02386 3279.0 + -dw 2.1e-9 + -Vm 1.39 28.3 0 -7.22 -0.59 # ref. 1 + Hnedkovsky et al., 1996, JCT 28, 125 +2H2S = (H2S)2 # activity correction for H2S solubility at high P, T + -analytical_expression 10.227 -0.01384 -2200 + -Vm 36.41 -71.95 0 0 2.58 +H2Sg = HSg- + H+ + -log_k -6.994 + -delta_h 5.30 kcal + -analytical_expression 11.17 -0.02386 -3279.0 + -gamma 3.5 0 + -dw 1.73e-9 + -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt +2H2Sg = (H2Sg)2 # activity correction for H2S solubility at high P, T + -analytical_expression 10.227 -0.01384 -2200 + -Vm 36.41 -71.95 0 0 2.58 +NO3- + 2 H+ + 2 e- = NO2- + H2O + -log_k 28.570 + -delta_h -43.760 kcal + -gamma 3.0 0 + -dw 1.91e-9 + -Vm 5.5864 5.8590 3.4472 -3.0212 1.1847 # supcrt +2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O + -log_k 207.08 + -delta_h -312.130 kcal + -dw 1.96e-9 + -Vm 7 # Pray et al., 1952, IEC 44. 1146 +NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O + -log_k 119.077 + -delta_h -187.055 kcal + -gamma 2.5 0 + -dw 1.98e-9 312 0.95 4.53 + -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # ref. 1 + +NH4+ = NH3 + H+ + -log_k -9.252 + -delta_h 12.48 kcal + -analytic 0.6322 -0.001225 -2835.76 + -dw 2.28e-9 + -Vm 6.69 2.8 3.58 -2.88 1.43 # ref. 2 +#NO3- + 10 H+ + 8 e- = AmmH+ + 3 H2O +# -log_k 119.077 +# -delta_h -187.055 kcal +# -gamma 2.5 0 +# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # ref. 1 + +#AmmH+ + SO4-2 = AmmHSO4- +NH4+ + SO4-2 = NH4SO4- + -log_k 1.11 + -Vm 14.0 0 -35.2 0 0 0 12.3 0 -0.141 1 # ref. 2 +H3BO3 = H2BO3- + H+ + -log_k -9.24 + -delta_h 3.224 kcal +H3BO3 + F- = BF(OH)3- + -log_k -0.4 + -delta_h 1.850 kcal +H3BO3 + 2 F- + H+ = BF2(OH)2- + H2O + -log_k 7.63 + -delta_h 1.618 kcal +H3BO3 + 2 H+ + 3 F- = BF3OH- + 2 H2O + -log_k 13.67 + -delta_h -1.614 kcal +H3BO3 + 3 H+ + 4 F- = BF4- + 3 H2O + -log_k 20.28 + -delta_h -1.846 kcal +PO4-3 + H+ = HPO4-2 + -log_k 12.346 + -delta_h -3.530 kcal + -gamma 5.0 0 + -dw 0.69e-9 + -Vm 3.52 1.09 8.39 -2.82 3.34 0 0 0 0 1 # ref. 2 +PO4-3 + 2 H+ = H2PO4- + -log_k 19.553 + -delta_h -4.520 kcal + -gamma 5.4 0 + -dw 0.846e-9 + -Vm 5.58 8.06 12.2 -3.11 1.3 0 0 0 1.62e-2 1 # ref. 2 +PO4-3 + 3H+ = H3PO4 + log_k 21.721 # log_k and delta_h from minteq.v4.dat, NIST46.3 + delta_h -10.1 kJ + -Vm 7.47 12.4 6.29 -3.29 0 # ref. 2 +H+ + F- = HF + -log_k 3.18 + -delta_h 3.18 kcal + -analytic -2.033 0.012645 429.01 + -Vm 3.4753 .7042 5.4732 -2.8081 -.0007 # supcrt +H+ + 2 F- = HF2- + -log_k 3.76 + -delta_h 4.550 kcal + -Vm 5.2263 4.9797 3.7928 -2.9849 1.2934 # supcrt +Ca+2 + H2O = CaOH+ + H+ + -log_k -12.78 +Ca+2 + CO3-2 = CaCO3 + -log_k 3.224 + -delta_h 3.545 kcal + -analytic -1228.732 -0.299440 35512.75 485.818 + -dw 4.46e-10 # complexes: calc'd with the Pikal formula + -Vm -.2430 -8.3748 9.0417 -2.4328 -.0300 # supcrt +Ca+2 + CO3-2 + H+ = CaHCO3+ + -log_k 11.435 + -delta_h -0.871 kcal + -analytic 1317.0071 0.34546894 -39916.84 -517.70761 563713.9 + -gamma 6.0 0 + -dw 5.06e-10 + -Vm 3.1911 .0104 5.7459 -2.7794 .3084 5.4 # supcrt +Ca+2 + SO4-2 = CaSO4 + -log_k 2.25 + -delta_h 1.325 kcal + -dw 4.71e-10 + -Vm 2.7910 -.9666 6.1300 -2.7390 -.0010 # supcrt +Ca+2 + HSO4- = CaHSO4+ + -log_k 1.08 +Ca+2 + PO4-3 = CaPO4- + -log_k 6.459 + -delta_h 3.10 kcal + -gamma 5.4 0.0 +Ca+2 + HPO4-2 = CaHPO4 + -log_k 2.739 + -delta_h 3.3 kcal +Ca+2 + H2PO4- = CaH2PO4+ + -log_k 1.408 + -delta_h 3.4 kcal + -gamma 5.4 0.0 +# Ca+2 + F- = CaF+ + # -log_k 0.94 + # -delta_h 4.120 kcal + # -gamma 5.5 0.0 + # -Vm .9846 -5.3773 7.8635 -2.5567 .6911 5.5 # supcrt +Mg+2 + H2O = MgOH+ + H+ + -log_k -11.44 + -delta_h 15.952 kcal + -gamma 6.5 0 +Mg+2 + CO3-2 = MgCO3 + -log_k 2.98 + -delta_h 2.713 kcal + -analytic 0.9910 0.00667 + -dw 4.21e-10 + -Vm -.5837 -9.2067 9.3687 -2.3984 -.0300 # supcrt +Mg+2 + H+ + CO3-2 = MgHCO3+ + -log_k 11.399 + -delta_h -2.771 kcal + -analytic 48.6721 0.03252849 -2614.335 -18.00263 563713.9 + -gamma 4.0 0 + -dw 4.78e-10 + -Vm 2.7171 -1.1469 6.2008 -2.7316 .5985 4 # supcrt +Mg+2 + SO4-2 = MgSO4 + -log_k 2.37 + -delta_h 4.550 kcal + -dw 4.45e-10 + -Vm 2.4 -0.97 6.1 -2.74 # est'd +Mg+2 + PO4-3 = MgPO4- + -log_k 6.589 + -delta_h 3.10 kcal + -gamma 5.4 0 +Mg+2 + HPO4-2 = MgHPO4 + -log_k 2.87 + -delta_h 3.3 kcal +Mg+2 + H2PO4- = MgH2PO4+ + -log_k 1.513 + -delta_h 3.4 kcal + -gamma 5.4 0 +Mg+2 + F- = MgF+ + -log_k 1.82 + -delta_h 3.20 kcal + -gamma 4.5 0 + -Vm .6494 -6.1958 8.1852 -2.5229 .9706 4.5 # supcrt +Na+ + OH- = NaOH + -log_k -10 # remove this complex +Na+ + CO3-2 = NaCO3- + -log_k 1.27 + -delta_h 8.91 kcal + -dw 1.2e-9 0 1e-10 1e-10 + -Vm 3.89 -8.23e-4 20 -9.44 3.02 9.05e-3 3.07 0 0.0233 1 # ref. 1 +Na+ + HCO3- = NaHCO3 + -log_k -0.25 + -delta_h -1 kcal + -dw 6.73e-10 + -Vm 0.431 # ref. 1 +Na+ + SO4-2 = NaSO4- + -log_k 0.7 + -delta_h 1.120 kcal + -gamma 5.4 0 + -dw 1.33e-9 0 0.57 1e-10 + -Vm 1e-5 16.4 -0.0678 -1.05 4.14 0 6.86 0 0.0242 0.53 # ref. 1 +Na+ + HPO4-2 = NaHPO4- + -log_k 0.29 + -gamma 5.4 0 + -Vm 5.2 8.1 13 -3 0.9 0 0 1.62e-2 1 # ref. 2 +Na+ + F- = NaF + -log_k -0.24 + -Vm 2.7483 -1.0708 6.1709 -2.7347 -.030 # supcrt +K+ + SO4-2 = KSO4- + -log_k 0.85 + -delta_h 2.250 kcal + -analytical 3.106 0.0 -673.6 + -gamma 5.4 0 + -dw 1.5e-9 0 1e-10 1e10 + -Vm 6.8 7.06 3.0 -2.07 1.1 0 0 0 0 1 # ref. 1 +K+ + HPO4-2 = KHPO4- + -log_k 0.29 + -gamma 5.4 0 + -Vm 5.4 8.1 19 -3.1 0.7 0 0 0 1.62e-2 1 # ref. 2 +Fe+2 + H2O = FeOH+ + H+ + -log_k -9.5 + -delta_h 13.20 kcal + -gamma 5.0 0 +Fe+2 + 3H2O = Fe(OH)3- + 3H+ + -log_k -31.0 + -delta_h 30.3 kcal + -gamma 5.0 0 +Fe+2 + Cl- = FeCl+ + -log_k 0.14 +Fe+2 + CO3-2 = FeCO3 + -log_k 4.38 +Fe+2 + HCO3- = FeHCO3+ + -log_k 2.0 +Fe+2 + SO4-2 = FeSO4 + -log_k 2.25 + -delta_h 3.230 kcal + -Vm -13 0 123 # ref. 2 +Fe+2 + HSO4- = FeHSO4+ + -log_k 1.08 +Fe+2 + 2HS- = Fe(HS)2 + -log_k 8.95 +Fe+2 + 3HS- = Fe(HS)3- + -log_k 10.987 +Fe+2 + HPO4-2 = FeHPO4 + -log_k 3.6 +Fe+2 + H2PO4- = FeH2PO4+ + -log_k 2.7 + -gamma 5.4 0 +Fe+2 + F- = FeF+ + -log_k 1.0 +Fe+2 = Fe+3 + e- + -log_k -13.02 + -delta_h 9.680 kcal + -gamma 9.0 0 +Fe+3 + H2O = FeOH+2 + H+ + -log_k -2.19 + -delta_h 10.4 kcal + -gamma 5.0 0 +Fe+3 + 2 H2O = Fe(OH)2+ + 2 H+ + -log_k -5.67 + -delta_h 17.1 kcal + -gamma 5.4 0 +Fe+3 + 3 H2O = Fe(OH)3 + 3 H+ + -log_k -12.56 + -delta_h 24.8 kcal +Fe+3 + 4 H2O = Fe(OH)4- + 4 H+ + -log_k -21.6 + -delta_h 31.9 kcal + -gamma 5.4 0 +Fe+2 + 2H2O = Fe(OH)2 + 2H+ + -log_k -20.57 + -delta_h 28.565 kcal +2 Fe+3 + 2 H2O = Fe2(OH)2+4 + 2 H+ + -log_k -2.95 + -delta_h 13.5 kcal +3 Fe+3 + 4 H2O = Fe3(OH)4+5 + 4 H+ + -log_k -6.3 + -delta_h 14.3 kcal +Fe+3 + Cl- = FeCl+2 + -log_k 1.48 + -delta_h 5.6 kcal + -gamma 5.0 0 +Fe+3 + 2 Cl- = FeCl2+ + -log_k 2.13 + -gamma 5.0 0 +Fe+3 + 3 Cl- = FeCl3 + -log_k 1.13 +Fe+3 + SO4-2 = FeSO4+ + -log_k 4.04 + -delta_h 3.91 kcal + -gamma 5.0 0 +Fe+3 + HSO4- = FeHSO4+2 + -log_k 2.48 +Fe+3 + 2 SO4-2 = Fe(SO4)2- + -log_k 5.38 + -delta_h 4.60 kcal +Fe+3 + HPO4-2 = FeHPO4+ + -log_k 5.43 + -delta_h 5.76 kcal + -gamma 5.0 0 +Fe+3 + H2PO4- = FeH2PO4+2 + -log_k 5.43 + -gamma 5.4 0 +Fe+3 + F- = FeF+2 + -log_k 6.2 + -delta_h 2.7 kcal + -gamma 5.0 0 +Fe+3 + 2 F- = FeF2+ + -log_k 10.8 + -delta_h 4.8 kcal + -gamma 5.0 0 +Fe+3 + 3 F- = FeF3 + -log_k 14.0 + -delta_h 5.4 kcal +Mn+2 + H2O = MnOH+ + H+ + -log_k -10.59 + -delta_h 14.40 kcal + -gamma 5.0 0 +Mn+2 + 3H2O = Mn(OH)3- + 3H+ + -log_k -34.8 + -gamma 5.0 0 +Mn+2 + Cl- = MnCl+ + -log_k 0.61 + -gamma 5.0 0 + -Vm 7.25 -1.08 -25.8 -2.73 3.99 5 0 0 0 1 # ref. 2 +Mn+2 + 2 Cl- = MnCl2 + -log_k 0.25 + -Vm 1e-5 0 144 # ref. 2 +Mn+2 + 3 Cl- = MnCl3- + -log_k -0.31 + -gamma 5.0 0 + -Vm 11.8 0 0 0 2.4 0 0 0 3.6e-2 1 # ref. 2 +Mn+2 + CO3-2 = MnCO3 + -log_k 4.9 +Mn+2 + HCO3- = MnHCO3+ + -log_k 1.95 + -gamma 5.0 0 +Mn+2 + SO4-2 = MnSO4 + -log_k 2.25 + -delta_h 3.370 kcal + -Vm -1.31 -1.83 62.3 -2.7 # ref. 2 +Mn+2 + 2 NO3- = Mn(NO3)2 + -log_k 0.6 + -delta_h -0.396 kcal + -Vm 6.16 0 29.4 0 0.9 # ref. 2 +Mn+2 + F- = MnF+ + -log_k 0.84 + -gamma 5.0 0 +Mn+2 = Mn+3 + e- + -log_k -25.51 + -delta_h 25.80 kcal + -gamma 9.0 0 +Al+3 + H2O = AlOH+2 + H+ + -log_k -5.0 + -delta_h 11.49 kcal + -analytic -38.253 0.0 -656.27 14.327 + -gamma 5.4 0 + -Vm -1.46 -11.4 10.2 -2.31 1.67 5.4 0 0 0 1 # ref. 2 and Barta and Hepler, 1986, Can. J. Chem. 64, 353. +Al+3 + 2 H2O = Al(OH)2+ + 2 H+ + -log_k -10.1 + -delta_h 26.90 kcal + -gamma 5.4 0 + -analytic 88.50 0.0 -9391.6 -27.121 +Al+3 + 3 H2O = Al(OH)3 + 3 H+ + -log_k -16.9 + -delta_h 39.89 kcal + -analytic 226.374 0.0 -18247.8 -73.597 +Al+3 + 4 H2O = Al(OH)4- + 4 H+ + -log_k -22.7 + -delta_h 42.30 kcal + -analytic 51.578 0.0 -11168.9 -14.865 + -gamma 4.5 0 + -dw 1.04e-9 # Mackin & Aller, 1983, GCA 47, 959 +Al+3 + SO4-2 = AlSO4+ + -log_k 3.5 + -delta_h 2.29 kcal + -gamma 4.5 0 +Al+3 + 2SO4-2 = Al(SO4)2- + -log_k 5.0 + -delta_h 3.11 kcal + -gamma 4.5 0 +Al+3 + HSO4- = AlHSO4+2 + -log_k 0.46 +Al+3 + F- = AlF+2 + -log_k 7.0 + -delta_h 1.060 kcal + -gamma 5.4 0 +Al+3 + 2 F- = AlF2+ + -log_k 12.7 + -delta_h 1.980 kcal + -gamma 5.4 0 +Al+3 + 3 F- = AlF3 + -log_k 16.8 + -delta_h 2.160 kcal +Al+3 + 4 F- = AlF4- + -log_k 19.4 + -delta_h 2.20 kcal + -gamma 4.5 0 +# Al+3 + 5 F- = AlF5-2 + # log_k 20.6 + # delta_h 1.840 kcal +# Al+3 + 6 F- = AlF6-3 + # log_k 20.6 + # delta_h -1.670 kcal +H4SiO4 = H3SiO4- + H+ + -log_k -9.83 + -delta_h 6.12 kcal + -analytic -302.3724 -0.050698 15669.69 108.18466 -1119669.0 + -gamma 4 0 + -Vm 7.94 1.0881 5.3224 -2.8240 1.4767 # supcrt + H2O in a1 +H4SiO4 = H2SiO4-2 + 2 H+ + -log_k -23.0 + -delta_h 17.6 kcal + -analytic -294.0184 -0.072650 11204.49 108.18466 -1119669.0 + -gamma 5.4 0 +H4SiO4 + 4 H+ + 6 F- = SiF6-2 + 4 H2O + -log_k 30.18 + -delta_h -16.260 kcal + -gamma 5.0 0 + -Vm 8.5311 13.0492 .6211 -3.3185 2.7716 # supcrt +Ba+2 + H2O = BaOH+ + H+ + -log_k -13.47 + -gamma 5.0 0 +Ba+2 + CO3-2 = BaCO3 + -log_k 2.71 + -delta_h 3.55 kcal + -analytic 0.113 0.008721 + -Vm .2907 -7.0717 8.5295 -2.4867 -.0300 # supcrt +Ba+2 + HCO3- = BaHCO3+ + -log_k 0.982 + -delta_h 5.56 kcal + -analytic -3.0938 0.013669 +Ba+2 + SO4-2 = BaSO4 + -log_k 2.7 +Sr+2 + H2O = SrOH+ + H+ + -log_k -13.29 + -gamma 5.0 0 +Sr+2 + CO3-2 + H+ = SrHCO3+ + -log_k 11.509 + -delta_h 2.489 kcal + -analytic 104.6391 0.04739549 -5151.79 -38.92561 563713.9 + -gamma 5.4 0 +Sr+2 + CO3-2 = SrCO3 + -log_k 2.81 + -delta_h 5.22 kcal + -analytic -1.019 0.012826 + -Vm -.1787 -8.2177 8.9799 -2.4393 -.0300 # supcrt +Sr+2 + SO4-2 = SrSO4 + -log_k 2.29 + -delta_h 2.08 kcal + -Vm 6.7910 -.9666 6.1300 -2.7390 -.0010 # celestite solubility +Li+ + SO4-2 = LiSO4- + -log_k 0.64 + -gamma 5.0 0 +Cu+2 + e- = Cu+ + -log_k 2.72 + -delta_h 1.65 kcal + -gamma 2.5 0 +Cu+ + 2Cl- = CuCl2- + -log_k 5.50 + -delta_h -0.42 kcal + -gamma 4.0 0 +Cu+ + 3Cl- = CuCl3-2 + -log_k 5.70 + -delta_h 0.26 kcal + -gamma 5.0 0.0 +Cu+2 + CO3-2 = CuCO3 + -log_k 6.73 +Cu+2 + 2CO3-2 = Cu(CO3)2-2 + -log_k 9.83 +Cu+2 + HCO3- = CuHCO3+ + -log_k 2.7 +Cu+2 + Cl- = CuCl+ + -log_k 0.43 + -delta_h 8.65 kcal + -gamma 4.0 0 + -Vm -4.19 0 30.4 0 0 4 0 0 1.94e-2 1 # ref. 2 +Cu+2 + 2Cl- = CuCl2 + -log_k 0.16 + -delta_h 10.56 kcal + -Vm 26.8 0 -136 # ref. 2 +Cu+2 + 3Cl- = CuCl3- + -log_k -2.29 + -delta_h 13.69 kcal + -gamma 4.0 0 +Cu+2 + 4Cl- = CuCl4-2 + -log_k -4.59 + -delta_h 17.78 kcal + -gamma 5.0 0 +Cu+2 + F- = CuF+ + -log_k 1.26 + -delta_h 1.62 kcal +Cu+2 + H2O = CuOH+ + H+ + -log_k -8.0 + -gamma 4.0 0 +Cu+2 + 2 H2O = Cu(OH)2 + 2 H+ + -log_k -13.68 +Cu+2 + 3 H2O = Cu(OH)3- + 3 H+ + -log_k -26.9 +Cu+2 + 4 H2O = Cu(OH)4-2 + 4 H+ + -log_k -39.6 +2Cu+2 + 2H2O = Cu2(OH)2+2 + 2H+ + -log_k -10.359 + -delta_h 17.539 kcal + -analytical 2.497 0.0 -3833.0 +Cu+2 + SO4-2 = CuSO4 + -log_k 2.31 + -delta_h 1.220 kcal + -Vm 5.21 0 -14.6 # ref. 2 +Cu+2 + 3HS- = Cu(HS)3- + -log_k 25.9 +Zn+2 + H2O = ZnOH+ + H+ + -log_k -8.96 + -delta_h 13.4 kcal +Zn+2 + 2 H2O = Zn(OH)2 + 2 H+ + -log_k -16.9 +Zn+2 + 3 H2O = Zn(OH)3- + 3 H+ + -log_k -28.4 +Zn+2 + 4 H2O = Zn(OH)4-2 + 4 H+ + -log_k -41.2 +Zn+2 + Cl- = ZnCl+ + -log_k 0.43 + -delta_h 7.79 kcal + -gamma 4.0 0 + -Vm 14.8 -3.91 -105.7 -2.62 0.203 4 0 0 -5.05e-2 1 # ref. 2 +Zn+2 + 2 Cl- = ZnCl2 + -log_k 0.45 + -delta_h 8.5 kcal + -Vm -10.1 4.57 241 -2.97 -1e-3 # ref. 2 +Zn+2 + 3Cl- = ZnCl3- + -log_k 0.5 + -delta_h 9.56 kcal + -gamma 4.0 0 + -Vm 0.772 15.5 -0.349 -3.42 1.25 0 -7.77 0 0 1 # ref. 2 +Zn+2 + 4Cl- = ZnCl4-2 + -log_k 0.2 + -delta_h 10.96 kcal + -gamma 5.0 0 + -Vm 28.42 28 -5.26 -3.94 2.67 0 0 0 4.62e-2 1 # ref. 2 +Zn+2 + H2O + Cl- = ZnOHCl + H+ + -log_k -7.48 +Zn+2 + 2HS- = Zn(HS)2 + -log_k 14.94 +Zn+2 + 3HS- = Zn(HS)3- + -log_k 16.1 +Zn+2 + CO3-2 = ZnCO3 + -log_k 5.3 +Zn+2 + 2CO3-2 = Zn(CO3)2-2 + -log_k 9.63 +Zn+2 + HCO3- = ZnHCO3+ + -log_k 2.1 +Zn+2 + SO4-2 = ZnSO4 + -log_k 2.37 + -delta_h 1.36 kcal + -Vm 2.51 0 18.8 # ref. 2 +Zn+2 + 2SO4-2 = Zn(SO4)2-2 + -log_k 3.28 + -Vm 10.9 0 -98.7 0 0 0 24 0 -0.236 1 # ref. 2 +Zn+2 + Br- = ZnBr+ + -log_k -0.58 +Zn+2 + 2Br- = ZnBr2 + -log_k -0.98 +Zn+2 + F- = ZnF+ + -log_k 1.15 + -delta_h 2.22 kcal +Cd+2 + H2O = CdOH+ + H+ + -log_k -10.08 + -delta_h 13.1 kcal +Cd+2 + 2 H2O = Cd(OH)2 + 2 H+ + -log_k -20.35 +Cd+2 + 3 H2O = Cd(OH)3- + 3 H+ + -log_k -33.3 +Cd+2 + 4 H2O = Cd(OH)4-2 + 4 H+ + -log_k -47.35 +2Cd+2 + H2O = Cd2OH+3 + H+ + -log_k -9.39 + -delta_h 10.9 kcal +Cd+2 + H2O + Cl- = CdOHCl + H+ + -log_k -7.404 + -delta_h 4.355 kcal +Cd+2 + NO3- = CdNO3+ + -log_k 0.4 + -delta_h -5.2 kcal + -Vm 5.95 0 -1.11 0 2.67 7 0 0 1.53e-2 1 # ref. 2 +Cd+2 + Cl- = CdCl+ + -log_k 1.98 + -delta_h 0.59 kcal + -Vm 5.69 0 -30.2 0 0 6 0 0 0.112 1 # ref. 2 +Cd+2 + 2 Cl- = CdCl2 + -log_k 2.6 + -delta_h 1.24 kcal + -Vm 5.53 # ref. 2 +Cd+2 + 3 Cl- = CdCl3- + -log_k 2.4 + -delta_h 3.9 kcal + -Vm 4.6 0 83.9 0 0 0 0 0 0 1 # ref. 2 +Cd+2 + CO3-2 = CdCO3 + -log_k 2.9 +Cd+2 + 2CO3-2 = Cd(CO3)2-2 + -log_k 6.4 +Cd+2 + HCO3- = CdHCO3+ + -log_k 1.5 +Cd+2 + SO4-2 = CdSO4 + -log_k 2.46 + -delta_h 1.08 kcal + -Vm 10.4 0 57.9 # ref. 2 +Cd+2 + 2SO4-2 = Cd(SO4)2-2 + -log_k 3.5 + -Vm -6.29 0 -93 0 9.5 7 0 0 0 1 # ref. 2 +Cd+2 + Br- = CdBr+ + -log_k 2.17 + -delta_h -0.81 kcal +Cd+2 + 2Br- = CdBr2 + -log_k 2.9 +Cd+2 + F- = CdF+ + -log_k 1.1 +Cd+2 + 2F- = CdF2 + -log_k 1.5 +Cd+2 + HS- = CdHS+ + -log_k 10.17 +Cd+2 + 2HS- = Cd(HS)2 + -log_k 16.53 +Cd+2 + 3HS- = Cd(HS)3- + -log_k 18.71 +Cd+2 + 4HS- = Cd(HS)4-2 + -log_k 20.9 +Pb+2 + H2O = PbOH+ + H+ + -log_k -7.71 +Pb+2 + 2 H2O = Pb(OH)2 + 2 H+ + -log_k -17.12 +Pb+2 + 3 H2O = Pb(OH)3- + 3 H+ + -log_k -28.06 +Pb+2 + 4 H2O = Pb(OH)4-2 + 4 H+ + -log_k -39.7 +2 Pb+2 + H2O = Pb2OH+3 + H+ + -log_k -6.36 +Pb+2 + Cl- = PbCl+ + -log_k 1.6 + -delta_h 4.38 kcal + -Vm 2.8934 -.7165 6.0316 -2.7494 .1281 6 # supcrt +Pb+2 + 2 Cl- = PbCl2 + -log_k 1.8 + -delta_h 1.08 kcal + -Vm 6.5402 8.1879 2.5318 -3.1175 -.0300 # supcrt +Pb+2 + 3 Cl- = PbCl3- + -log_k 1.7 + -delta_h 2.17 kcal + -Vm 11.0396 19.1743 -1.7863 -3.5717 .7356 # supcrt +Pb+2 + 4 Cl- = PbCl4-2 + -log_k 1.38 + -delta_h 3.53 kcal + -Vm 16.4150 32.2997 -6.9452 -4.1143 2.3118 # supcrt +Pb+2 + CO3-2 = PbCO3 + -log_k 7.24 +Pb+2 + 2 CO3-2 = Pb(CO3)2-2 + -log_k 10.64 +Pb+2 + HCO3- = PbHCO3+ + -log_k 2.9 +Pb+2 + SO4-2 = PbSO4 + -log_k 2.75 +Pb+2 + 2 SO4-2 = Pb(SO4)2-2 + -log_k 3.47 +Pb+2 + 2HS- = Pb(HS)2 + -log_k 15.27 +Pb+2 + 3HS- = Pb(HS)3- + -log_k 16.57 +3Pb+2 + 4H2O = Pb3(OH)4+2 + 4H+ + -log_k -23.88 + -delta_h 26.5 kcal +Pb+2 + NO3- = PbNO3+ + -log_k 1.17 +Pb+2 + Br- = PbBr+ + -log_k 1.77 + -delta_h 2.88 kcal +Pb+2 + 2Br- = PbBr2 + -log_k 1.44 +Pb+2 + F- = PbF+ + -log_k 1.25 +Pb+2 + 2F- = PbF2 + -log_k 2.56 +Pb+2 + 3F- = PbF3- + -log_k 3.42 +Pb+2 + 4F- = PbF4-2 + -log_k 3.1 + +PHASES +Calcite + CaCO3 = CO3-2 + Ca+2 + -log_k -8.48 + -delta_h -2.297 kcal + -analytic 17.118 -0.046528 -3496 # 0 - 250°C, Ellis, 1959, Plummer and Busenberg, 1982 + -Vm 36.9 cm3/mol # MW (100.09 g/mol) / rho (2.71 g/cm3) +Aragonite + CaCO3 = CO3-2 + Ca+2 + -log_k -8.336 + -delta_h -2.589 kcal + -analytic -171.9773 -0.077993 2903.293 71.595 + -Vm 34.04 +Dolomite + CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2 + -log_k -17.09 + -delta_h -9.436 kcal + -analytic 31.283 -0.0898 -6438 # 25°C: Hemingway and Robie, 1994; 50–175°C: Bénézeth et al., 2018, GCA 224, 262-275. + -Vm 64.5 +Siderite + FeCO3 = Fe+2 + CO3-2 + -log_k -10.89 + -delta_h -2.480 kcal + -Vm 29.2 +Rhodochrosite + MnCO3 = Mn+2 + CO3-2 + -log_k -11.13 + -delta_h -1.430 kcal + -Vm 31.1 +Strontianite + SrCO3 = Sr+2 + CO3-2 + -log_k -9.271 + -delta_h -0.400 kcal + -analytic 155.0305 0.0 -7239.594 -56.58638 + -Vm 39.69 +Witherite + BaCO3 = Ba+2 + CO3-2 + -log_k -8.562 + -delta_h 0.703 kcal + -analytic 607.642 0.121098 -20011.25 -236.4948 + -Vm 46 +Gypsum + CaSO4:2H2O = Ca+2 + SO4-2 + 2 H2O + -log_k -4.58 + -delta_h -0.109 kcal + -analytic 68.2401 0.0 -3221.51 -25.0627 + -analytical_expression 93.7 5.99E-03 -4e3 -35.019 # better fits the appendix data of Appelo, 2015, AG 55, 62 + -Vm 73.9 # 172.18 / 2.33 (Vm H2O = 13.9 cm3/mol) +Anhydrite + CaSO4 = Ca+2 + SO4-2 + -log_k -4.36 + -delta_h -1.710 kcal + -analytic 84.90 0 -3135.12 -31.79 # 50 - 160oC, 1 - 1e3 atm, anhydrite dissolution, Blount and Dickson, 1973, Am. Mineral. 58, 323. + -Vm 46.1 # 136.14 / 2.95 +Celestite + SrSO4 = Sr+2 + SO4-2 + -log_k -6.63 + -delta_h -4.037 kcal +# -analytic -14805.9622 -2.4660924 756968.533 5436.3588 -40553604.0 + -analytic -7.14 6.11e-3 75 0 0 -1.79e-5 # Howell et al., 1992, JCED 37, 464. + -Vm 46.4 +Barite + BaSO4 = Ba+2 + SO4-2 + -log_k -9.97 + -delta_h 6.35 kcal + -analytical_expression -282.43 -8.972e-2 5822 113.08 # Blount 1977; Templeton, 1960 + -Vm 52.9 +Hydroxyapatite + Ca5(PO4)3OH + 4 H+ = H2O + 3 HPO4-2 + 5 Ca+2 + -log_k -3.421 + -delta_h -36.155 kcal + -Vm 128.9 +Fluorite + CaF2 = Ca+2 + 2 F- + -log_k -10.6 + -delta_h 4.69 kcal + -analytic 66.348 0.0 -4298.2 -25.271 + -Vm 15.7 +SiO2(a) + SiO2 + 2 H2O = H4SiO4 + -log_k -2.71 + -delta_h 3.340 kcal + -analytic -0.26 0.0 -731.0 +Chalcedony + SiO2 + 2 H2O = H4SiO4 + -log_k -3.55 + -delta_h 4.720 kcal + -analytic -0.09 0.0 -1032.0 + -Vm 23.1 +Quartz + SiO2 + 2 H2O = H4SiO4 + -log_k -3.98 + -delta_h 5.990 kcal + -analytic 0.41 0.0 -1309.0 + -Vm 22.67 +Gibbsite + Al(OH)3 + 3 H+ = Al+3 + 3 H2O + -log_k 8.11 + -delta_h -22.800 kcal + -Vm 32.22 +Al(OH)3(a) + Al(OH)3 + 3 H+ = Al+3 + 3 H2O + -log_k 10.8 + -delta_h -26.500 kcal +Kaolinite + Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3 + -log_k 7.435 + -delta_h -35.300 kcal + -Vm 99.35 +Albite + NaAlSi3O8 + 8 H2O = Na+ + Al(OH)4- + 3 H4SiO4 + -log_k -18.002 + -delta_h 25.896 kcal + -Vm 101.31 +Anorthite + CaAl2Si2O8 + 8 H2O = Ca+2 + 2 Al(OH)4- + 2 H4SiO4 + -log_k -19.714 + -delta_h 11.580 kcal + -Vm 105.05 +K-feldspar + KAlSi3O8 + 8 H2O = K+ + Al(OH)4- + 3 H4SiO4 + -log_k -20.573 + -delta_h 30.820 kcal + -Vm 108.15 +K-mica + KAl3Si3O10(OH)2 + 10 H+ = K+ + 3 Al+3 + 3 H4SiO4 + -log_k 12.703 + -delta_h -59.376 kcal +Chlorite(14A) + Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3H4SiO4 + 6H2O + -log_k 68.38 + -delta_h -151.494 kcal +Ca-Montmorillonite + Ca0.165Al2.33Si3.67O10(OH)2 + 12 H2O = 0.165Ca+2 + 2.33 Al(OH)4- + 3.67 H4SiO4 + 2 H+ + -log_k -45.027 + -delta_h 58.373 kcal + -Vm 156.16 +Talc + Mg3Si4O10(OH)2 + 4 H2O + 6 H+ = 3 Mg+2 + 4 H4SiO4 + -log_k 21.399 + -delta_h -46.352 kcal + -Vm 68.34 +Illite + K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 0.6K+ + 0.25Mg+2 + 2.3Al(OH)4- + 3.5H4SiO4 + 1.2H+ + -log_k -40.267 + -delta_h 54.684 kcal + -Vm 141.48 +Chrysotile + Mg3Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 3 Mg+2 + -log_k 32.2 + -delta_h -46.800 kcal + -analytic 13.248 0.0 10217.1 -6.1894 + -Vm 106.5808 # 277.11/2.60 +Sepiolite + Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 + -log_k 15.760 + -delta_h -10.700 kcal + -Vm 143.765 +Sepiolite(d) + Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 + -log_k 18.66 +Hematite + Fe2O3 + 6 H+ = 2 Fe+3 + 3 H2O + -log_k -4.008 + -delta_h -30.845 kcal + -Vm 30.39 +Goethite + FeOOH + 3 H+ = Fe+3 + 2 H2O + -log_k -1.0 + -delta_h -14.48 kcal + -Vm 20.84 +Fe(OH)3(a) + Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O + -log_k 4.891 +Pyrite + FeS2 + 2 H+ + 2 e- = Fe+2 + 2 HS- + -log_k -18.479 + -delta_h 11.300 kcal + -Vm 23.48 +FeS(ppt) + FeS + H+ = Fe+2 + HS- + -log_k -3.915 +Mackinawite + FeS + H+ = Fe+2 + HS- + -log_k -4.648 + -Vm 20.45 +Sulfur + S + 2H+ + 2e- = H2S + -log_k 4.882 + -delta_h -9.5 kcal +Vivianite + Fe3(PO4)2:8H2O = 3 Fe+2 + 2 PO4-3 + 8 H2O + -log_k -36.0 +Pyrolusite # H2O added for surface calc's + MnO2:H2O + 4 H+ + 2 e- = Mn+2 + 3 H2O + -log_k 41.38 + -delta_h -65.110 kcal +Hausmannite + Mn3O4 + 8 H+ + 2 e- = 3 Mn+2 + 4 H2O + -log_k 61.03 + -delta_h -100.640 kcal +Manganite + MnOOH + 3 H+ + e- = Mn+2 + 2 H2O + -log_k 25.34 +Pyrochroite + Mn(OH)2 + 2 H+ = Mn+2 + 2 H2O + -log_k 15.2 +Halite + NaCl = Cl- + Na+ + log_k 1.570 + -delta_h 1.37 + #-analytic -713.4616 -.1201241 37302.21 262.4583 -2106915. + -Vm 27.1 +Sylvite + KCl = K+ + Cl- + log_k 0.900 + -delta_h 8.5 + # -analytic 3.984 0.0 -919.55 + Vm 37.5 +# Gases... +CO2(g) + CO2 = CO2 + -log_k -1.468 + -delta_h -4.776 kcal + -analytic 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5 + -T_c 304.2 # critical T, K + -P_c 72.86 # critical P, atm + -Omega 0.225 # acentric factor +H2O(g) + H2O = H2O + -log_k 1.506; delta_h -44.03 kJ + -T_c 647.3 + -P_c 217.60 + -Omega 0.344 + -analytic -16.5066 -2.0013E-3 2710.7 3.7646 0 2.24E-6 +O2(g) + O2 = O2 + -log_k -2.8983 + -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 + -T_c 154.6; -P_c 49.80; -Omega 0.021 +H2(g) + H2 = H2 + -log_k -3.1050 + -delta_h -4.184 kJ + -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 + -T_c 33.2; -P_c 12.80; -Omega -0.225 +N2(g) + N2 = N2 + -log_k -3.1864 + -analytic -58.453 1.818e-3 3199 17.909 -27460 + -T_c 126.2; -P_c 33.50; -Omega 0.039 +H2S(g) + H2S = H+ + HS- + log_k -7.93 + -delta_h 9.1 + -analytic -45.07 -0.02418 0 17.9205 # H2S solubilities, 0 - 300°C, 1 - 987 atm, Jiang et al., 2020, CG 555, 119816 + -T_c 373.2; -P_c 88.20; -Omega 0.1 +CH4(g) + CH4 = CH4 + -log_k -2.8 + -analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100°C + -T_c 190.6 ; -P_c 45.40 ; -Omega 0.008 +#Amm(g) +# Amm = Amm +NH3(g) + NH3 = NH3 + -log_k 1.7966 + -analytic -18.758 3.3670e-4 2.5113e3 4.8619 39.192 + -T_c 405.6; -P_c 111.3; -Omega 0.25 +# redox-uncoupled gases +Oxg(g) + Oxg = Oxg + -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 + -T_c 154.6 ; -P_c 49.80 ; -Omega 0.021 +Hdg(g) + Hdg = Hdg + -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 + -T_c 33.2 ; -P_c 12.80 ; -Omega -0.225 +Ntg(g) + Ntg = Ntg + -analytic -58.453 1.81800e-3 3199 17.909 -27460 + T_c 126.2 ; -P_c 33.50 ; -Omega 0.039 +Mtg(g) + Mtg = Mtg + -log_k -2.8 + -analytic 10.44 -7.65e-3 -6669 0 1.014e6 # CH4 solubilities 25 - 100°C + -T_c 190.6 ; -P_c 45.40 ; -Omega 0.008 +H2Sg(g) + H2Sg = H+ + HSg- + log_k -7.93 + -delta_h 9.1 + -analytic -45.07 -0.02418 0 17.9205 # H2S solubilities, 0 - 300°C, 1 - 987 atm, Jiang et al., 2020, CG 555, 119816 + -T_c 373.2 ; -P_c 88.20 ; -Omega 0.1 +Melanterite + FeSO4:7H2O = 7 H2O + Fe+2 + SO4-2 + -log_k -2.209 + -delta_h 4.910 kcal + -analytic 1.447 -0.004153 0.0 0.0 -214949.0 +Alunite + KAl3(SO4)2(OH)6 + 6 H+ = K+ + 3 Al+3 + 2 SO4-2 + 6H2O + -log_k -1.4 + -delta_h -50.250 kcal +Jarosite-K + KFe3(SO4)2(OH)6 + 6 H+ = 3 Fe+3 + 6 H2O + K+ + 2 SO4-2 + -log_k -9.21 + -delta_h -31.280 kcal +Zn(OH)2(e) + Zn(OH)2 + 2 H+ = Zn+2 + 2 H2O + -log_k 11.5 +Smithsonite + ZnCO3 = Zn+2 + CO3-2 + -log_k -10.0 + -delta_h -4.36 kcal +Sphalerite + ZnS + H+ = Zn+2 + HS- + -log_k -11.618 + -delta_h 8.250 kcal +Willemite 289 + Zn2SiO4 + 4H+ = 2Zn+2 + H4SiO4 + -log_k 15.33 + -delta_h -33.37 kcal +Cd(OH)2 + Cd(OH)2 + 2 H+ = Cd+2 + 2 H2O + -log_k 13.65 +Otavite 315 + CdCO3 = Cd+2 + CO3-2 + -log_k -12.1 + -delta_h -0.019 kcal +CdSiO3 328 + CdSiO3 + H2O + 2H+ = Cd+2 + H4SiO4 + -log_k 9.06 + -delta_h -16.63 kcal +CdSO4 329 + CdSO4 = Cd+2 + SO4-2 + -log_k -0.1 + -delta_h -14.74 kcal +Cerussite 365 + PbCO3 = Pb+2 + CO3-2 + -log_k -13.13 + -delta_h 4.86 kcal +Anglesite 384 + PbSO4 = Pb+2 + SO4-2 + -log_k -7.79 + -delta_h 2.15 kcal +Pb(OH)2 389 + Pb(OH)2 + 2H+ = Pb+2 + 2H2O + -log_k 8.15 + -delta_h -13.99 kcal + +EXCHANGE_MASTER_SPECIES + X X- +EXCHANGE_SPECIES + X- = X- + -log_k 0.0 + + Na+ + X- = NaX + -log_k 0.0 + -gamma 4.08 0.082 + + K+ + X- = KX + -log_k 0.7 + -gamma 3.5 0.015 + -delta_h -4.3 # Jardine & Sparks, 1984 + + Li+ + X- = LiX + -log_k -0.08 + -gamma 6.0 0 + -delta_h 1.4 # Merriam & Thomas, 1956 + +# !!!!! +# H+ + X- = HX +# -log_k 1.0 +# -gamma 9.0 0 + +# AmmH+ + X- = AmmHX + NH4+ + X- = NH4X + -log_k 0.6 + -gamma 2.5 0 + -delta_h -2.4 # Laudelout et al., 1968 + + Ca+2 + 2X- = CaX2 + -log_k 0.8 + -gamma 5.0 0.165 + -delta_h 7.2 # Van Bladel & Gheyl, 1980 + + Mg+2 + 2X- = MgX2 + -log_k 0.6 + -gamma 5.5 0.2 + -delta_h 7.4 # Laudelout et al., 1968 + + Sr+2 + 2X- = SrX2 + -log_k 0.91 + -gamma 5.26 0.121 + -delta_h 5.5 # Laudelout et al., 1968 + + Ba+2 + 2X- = BaX2 + -log_k 0.91 + -gamma 4.0 0.153 + -delta_h 4.5 # Laudelout et al., 1968 + + Mn+2 + 2X- = MnX2 + -log_k 0.52 + -gamma 6.0 0 + + Fe+2 + 2X- = FeX2 + -log_k 0.44 + -gamma 6.0 0 + + Cu+2 + 2X- = CuX2 + -log_k 0.6 + -gamma 6.0 0 + + Zn+2 + 2X- = ZnX2 + -log_k 0.8 + -gamma 5.0 0 + + Cd+2 + 2X- = CdX2 + -log_k 0.8 + -gamma 0.0 0 + + Pb+2 + 2X- = PbX2 + -log_k 1.05 + -gamma 0.0 0 + + Al+3 + 3X- = AlX3 + -log_k 0.41 + -gamma 9.0 0 + + AlOH+2 + 2X- = AlOHX2 + -log_k 0.89 + -gamma 0.0 0 + +SURFACE_MASTER_SPECIES + Hfo_s Hfo_sOH + Hfo_w Hfo_wOH +SURFACE_SPECIES +# All surface data from +# Dzombak and Morel, 1990 +# +# +# Acid-base data from table 5.7 +# +# strong binding site--Hfo_s, + + Hfo_sOH = Hfo_sOH + -log_k 0 + + Hfo_sOH + H+ = Hfo_sOH2+ + -log_k 7.29 # = pKa1,int + + Hfo_sOH = Hfo_sO- + H+ + -log_k -8.93 # = -pKa2,int + +# weak binding site--Hfo_w + + Hfo_wOH = Hfo_wOH + -log_k 0 + + Hfo_wOH + H+ = Hfo_wOH2+ + -log_k 7.29 # = pKa1,int + + Hfo_wOH = Hfo_wO- + H+ + -log_k -8.93 # = -pKa2,int +############################################### +# CATIONS # +############################################### +# +# Cations from table 10.1 or 10.5 +# +# Calcium + Hfo_sOH + Ca+2 = Hfo_sOHCa+2 + -log_k 4.97 + + Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+ + -log_k -5.85 +# Strontium + Hfo_sOH + Sr+2 = Hfo_sOHSr+2 + -log_k 5.01 + + Hfo_wOH + Sr+2 = Hfo_wOSr+ + H+ + -log_k -6.58 + + Hfo_wOH + Sr+2 + H2O = Hfo_wOSrOH + 2H+ + -log_k -17.6 +# Barium + Hfo_sOH + Ba+2 = Hfo_sOHBa+2 + -log_k 5.46 + + Hfo_wOH + Ba+2 = Hfo_wOBa+ + H+ + -log_k -7.2 # table 10.5 +# +# Cations from table 10.2 +# +# Cadmium + Hfo_sOH + Cd+2 = Hfo_sOCd+ + H+ + -log_k 0.47 + + Hfo_wOH + Cd+2 = Hfo_wOCd+ + H+ + -log_k -2.91 +# Zinc + Hfo_sOH + Zn+2 = Hfo_sOZn+ + H+ + -log_k 0.99 + + Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+ + -log_k -1.99 +# Copper + Hfo_sOH + Cu+2 = Hfo_sOCu+ + H+ + -log_k 2.89 + + Hfo_wOH + Cu+2 = Hfo_wOCu+ + H+ + -log_k 0.6 # table 10.5 +# Lead + Hfo_sOH + Pb+2 = Hfo_sOPb+ + H+ + -log_k 4.65 + + Hfo_wOH + Pb+2 = Hfo_wOPb+ + H+ + -log_k 0.3 # table 10.5 +# +# Derived constants table 10.5 +# +# Magnesium + Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+ + -log_k -4.6 +# Manganese + Hfo_sOH + Mn+2 = Hfo_sOMn+ + H+ + -log_k -0.4 # table 10.5 + + Hfo_wOH + Mn+2 = Hfo_wOMn+ + H+ + -log_k -3.5 # table 10.5 +# Iron, strong site: Appelo, Van der Weiden, Tournassat & Charlet, EST 36, 3096 + Hfo_sOH + Fe+2 = Hfo_sOFe+ + H+ + -log_k -0.95 +# Iron, weak site: Liger et al., GCA 63, 2939, re-optimized for D&M + Hfo_wOH + Fe+2 = Hfo_wOFe+ + H+ + -log_k -2.98 + + Hfo_wOH + Fe+2 + H2O = Hfo_wOFeOH + 2H+ + -log_k -11.55 +############################################### +# ANIONS # +############################################### +# +# Anions from table 10.6 +# +# Phosphate + Hfo_wOH + PO4-3 + 3H+ = Hfo_wH2PO4 + H2O + -log_k 31.29 + + Hfo_wOH + PO4-3 + 2H+ = Hfo_wHPO4- + H2O + -log_k 25.39 + + Hfo_wOH + PO4-3 + H+ = Hfo_wPO4-2 + H2O + -log_k 17.72 +# +# Anions from table 10.7 +# +# Borate + Hfo_wOH + H3BO3 = Hfo_wH2BO3 + H2O + -log_k 0.62 +# +# Anions from table 10.8 +# +# Sulfate + Hfo_wOH + SO4-2 + H+ = Hfo_wSO4- + H2O + -log_k 7.78 + + Hfo_wOH + SO4-2 = Hfo_wOHSO4-2 + -log_k 0.79 +# +# Derived constants table 10.10 +# + Hfo_wOH + F- + H+ = Hfo_wF + H2O + -log_k 8.7 + + Hfo_wOH + F- = Hfo_wOHF- + -log_k 1.6 +# +# Carbonate: Van Geen et al., 1994 reoptimized for D&M model +# + Hfo_wOH + CO3-2 + H+ = Hfo_wCO3- + H2O + -log_k 12.56 + + Hfo_wOH + CO3-2 + 2H+= Hfo_wHCO3 + H2O + -log_k 20.62 +# +# Silicate: Swedlund, P.J. and Webster, J.G., 1999. Water Research 33, 3413-3422. +# + Hfo_wOH + H4SiO4 = Hfo_wH3SiO4 + H2O ; log_K 4.28 + Hfo_wOH + H4SiO4 = Hfo_wH2SiO4- + H+ + H2O ; log_K -3.22 + Hfo_wOH + H4SiO4 = Hfo_wHSiO4-2 + 2H+ + H2O ; log_K -11.69 + +RATES + +########### +#Quartz +########### +# +####### +# Example of quartz kinetic rates block: +# KINETICS +# Quartz +# -m0 158.8 # 90 % Qu +# -parms 0.146 1.5 +# -step 3.1536e8 in 10 +# -tol 1e-12 + +Quartz + -start +1 REM Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683 +2 REM k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol +3 REM sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259) +4 REM PARM(1) = Specific area of Quartz, m^2/mol Quartz +5 REM PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35 + +10 dif_temp = 1/TK - 1/298 +20 pk_w = 13.7 + 4700.4 * dif_temp +40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 - SR("Quartz")) +# Integrate... +50 SAVE moles * TIME + -end + +########### +#K-feldspar +########### +# +# Sverdrup and Warfvinge, 1995, Estimating field weathering rates +# using laboratory kinetics: Reviews in mineralogy and geochemistry, +# vol. 31, p. 485-541. +# +# As described in: +# Appelo and Postma, 2005, Geochemistry, groundwater +# and pollution, 2nd Edition: A.A. Balkema Publishers, +# p. 162-163 and 395-399. +# +# Assume soil is 10% K-feldspar by mass in 1 mm spheres (radius 0.05 mm) +# Assume density of rock and Kspar is 2600 kg/m^3 = 2.6 kg/L +# GFW Kspar 0.278 kg/mol +# +# Moles of Kspar per liter pore space calculation: +# Mass of rock per liter pore space = 0.7*2.6/0.3 = 6.07 kg rock/L pore space +# Mass of Kspar per liter pore space 6.07x0.1 = 0.607 kg Kspar/L pore space +# Moles of Kspar per liter pore space 0.607/0.278 = 2.18 mol Kspar/L pore space +# +# Specific area calculation: +# Volume of sphere 4/3 x pi x r^3 = 5.24e-13 m^3 Kspar/sphere +# Mass of sphere 2600 x 5.24e-13 = 1.36e-9 kg Kspar/sphere +# Moles of Kspar in sphere 1.36e-9/0.278 = 4.90e-9 mol Kspar/sphere +# Surface area of one sphere 4 x pi x r^2 = 3.14e-8 m^2/sphere +# Specific area of K-feldspar in sphere 3.14e-8/4.90e-9 = 6.41 m^2/mol Kspar +# +# +# Example of KINETICS data block for K-feldspar rate: +# KINETICS 1 +# K-feldspar +# -m0 2.18 # 10% Kspar, 0.1 mm cubes +# -m 2.18 # Moles per L pore space +# -parms 6.41 0.1 # m^2/mol Kspar, fraction adjusts lab rate to field rate +# -time 1.5 year in 40 + +K-feldspar + -start +1 REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1 +2 REM PARM(1) = Specific area of Kspar m^2/mol Kspar +3 REM PARM(2) = Adjusts lab rate to field rate +4 REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281) +5 REM K-Feldspar parameters +10 DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3 +20 RESTORE 10 +30 READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH +40 DATA 3500, 2000, 2500, 2000 +50 RESTORE 40 +60 READ e_H, e_H2O, e_OH, e_CO2 +70 pk_CO2 = 13 +80 n_CO2 = 0.6 +100 REM Generic rate follows +110 dif_temp = 1/TK - 1/281 +120 BC = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2") +130 REM rate by H+ +140 pk_H = pk_H + e_H * dif_temp +150 rate_H = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC) +160 REM rate by hydrolysis +170 pk_H2O = pk_H2O + e_H2O * dif_temp +180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC) +190 REM rate by OH- +200 pk_OH = pk_OH + e_OH * dif_temp +210 rate_OH = 10^-pk_OH * ACT("OH-")^o_OH +220 REM rate by CO2 +230 pk_CO2 = pk_CO2 + e_CO2 * dif_temp +240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2 +250 rate = rate_H + rate_H2O + rate_OH + rate_CO2 +260 area = PARM(1) * M0 *(M/M0)^0.67 +270 rate = PARM(2) * area * rate * (1-SR("K-feldspar")) +280 moles = rate * TIME +290 SAVE moles + -end + + +########### +#Albite +########### +# +# Sverdrup and Warfvinge, 1995, Estimating field weathering rates +# using laboratory kinetics: Reviews in mineralogy and geochemistry, +# vol. 31, p. 485-541. +# +# As described in: +# Appelo and Postma, 2005, Geochemistry, groundwater +# and pollution, 2nd Edition: A.A. Balkema Publishers, +# p. 162-163 and 395-399. +# +# Example of KINETICS data block for Albite rate: +# KINETICS 1 +# Albite +# -m0 0.46 # 2% Albite, 0.1 mm cubes +# -m 0.46 # Moles per L pore space +# -parms 6.04 0.1 # m^2/mol Albite, fraction adjusts lab rate to field rate +# -time 1.5 year in 40 +# +# Assume soil is 2% Albite by mass in 1 mm spheres (radius 0.05 mm) +# Assume density of rock and Albite is 2600 kg/m^3 = 2.6 kg/L +# GFW Albite 0.262 kg/mol +# +# Moles of Albite per liter pore space calculation: +# Mass of rock per liter pore space = 0.7*2.6/0.3 = 6.07 kg rock/L pore space +# Mass of Albite per liter pore space 6.07x0.02 = 0.121 kg Albite/L pore space +# Moles of Albite per liter pore space 0.607/0.262 = 0.46 mol Albite/L pore space +# +# Specific area calculation: +# Volume of sphere 4/3 x pi x r^3 = 5.24e-13 m^3 Albite/sphere +# Mass of sphere 2600 x 5.24e-13 = 1.36e-9 kg Albite/sphere +# Moles of Albite in sphere 1.36e-9/0.262 = 5.20e-9 mol Albite/sphere +# Surface area of one sphere 4 x pi x r^2 = 3.14e-8 m^2/sphere +# Specific area of Albite in sphere 3.14e-8/5.20e-9 = 6.04 m^2/mol Albite + +Albite + -start +1 REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1 +2 REM PARM(1) = Specific area of Albite m^2/mol Albite +3 REM PARM(2) = Adjusts lab rate to field rate +4 REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281) +5 REM Albite parameters +10 DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3 +20 RESTORE 10 +30 READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH +40 DATA 3500, 2000, 2500, 2000 +50 RESTORE 40 +60 READ e_H, e_H2O, e_OH, e_CO2 +70 pk_CO2 = 13 +80 n_CO2 = 0.6 +100 REM Generic rate follows +110 dif_temp = 1/TK - 1/281 +120 BC = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2") +130 REM rate by H+ +140 pk_H = pk_H + e_H * dif_temp +150 rate_H = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC) +160 REM rate by hydrolysis +170 pk_H2O = pk_H2O + e_H2O * dif_temp +180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC) +190 REM rate by OH- +200 pk_OH = pk_OH + e_OH * dif_temp +210 rate_OH = 10^-pk_OH * ACT("OH-")^o_OH +220 REM rate by CO2 +230 pk_CO2 = pk_CO2 + e_CO2 * dif_temp +240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2 +250 rate = rate_H + rate_H2O + rate_OH + rate_CO2 +260 area = PARM(1) * M0 *(M/M0)^0.67 +270 rate = PARM(2) * area * rate * (1-SR("Albite")) +280 moles = rate * TIME +290 SAVE moles + -end + +######## +#Calcite +######## +# Example of KINETICS data block for calcite rate, +# in mmol/cm2/s, Plummer et al., 1978, AJS 278, 179; Appelo et al., AG 13, 257. +# KINETICS 1 +# Calcite +# -tol 1e-8 +# -m0 3.e-3 +# -m 3.e-3 +# -parms 1.67e5 0.6 # cm^2/mol calcite, exp factor +# -time 1 day + +Calcite + -start +1 REM PARM(1) = specific surface area of calcite, cm^2/mol calcite +2 REM PARM(2) = exponent for M/M0 + +10 si_cc = SI("Calcite") +20 IF (M <= 0 and si_cc < 0) THEN GOTO 200 +30 k1 = 10^(0.198 - 444.0 / TK ) +40 k2 = 10^(2.84 - 2177.0 /TK ) +50 IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK) +60 IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK ) +80 IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M +110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O")) +120 rate = rate * (1 - 10^(2/3*si_cc)) +130 moles = rate * 0.001 * TIME # convert from mmol to mol +200 SAVE moles + -end + +####### +#Pyrite +####### +# +# Williamson, M.A. and Rimstidt, J.D., 1994, +# Geochimica et Cosmochimica Acta, v. 58, p. 5443-5454, +# rate equation is mol m^-2 s^-1. +# +# Example of KINETICS data block for pyrite rate: +# KINETICS 1 +# Pyrite +# -tol 1e-8 +# -m0 5.e-4 +# -m 5.e-4 +# -parms 0.3 0.67 .5 -0.11 +# -time 1 day in 10 +Pyrite + -start +1 REM Williamson and Rimstidt, 1994 +2 REM PARM(1) = log10(specific area), log10(m^2 per mole pyrite) +3 REM PARM(2) = exp for (M/M0) +4 REM PARM(3) = exp for O2 +5 REM PARM(4) = exp for H+ + +10 REM Dissolution in presence of DO +20 if (M <= 0) THEN GOTO 200 +30 if (SI("Pyrite") >= 0) THEN GOTO 200 +40 log_rate = -8.19 + PARM(3)*LM("O2") + PARM(4)*LM("H+") +50 log_area = PARM(1) + LOG10(M0) + PARM(2)*LOG10(M/M0) +60 moles = 10^(log_area + log_rate) * TIME +200 SAVE moles + -end + +########## +#Organic_C +########## +# +# Example of KINETICS data block for SOC (sediment organic carbon): +# KINETICS 1 +# Organic_C +# -formula C +# -tol 1e-8 +# -m 5e-3 # SOC in mol +# -time 30 year in 15 +Organic_C + -start +1 REM Additive Monod kinetics for SOC (sediment organic carbon) +2 REM Electron acceptors: O2, NO3, and SO4 + +10 if (M <= 0) THEN GOTO 200 +20 mO2 = MOL("O2") +30 mNO3 = TOT("N(5)") +40 mSO4 = TOT("S(6)") +50 k_O2 = 1.57e-9 # 1/sec +60 k_NO3 = 1.67e-11 # 1/sec +70 k_SO4 = 1.e-13 # 1/sec +80 rate = k_O2 * mO2/(2.94e-4 + mO2) +90 rate = rate + k_NO3 * mNO3/(1.55e-4 + mNO3) +100 rate = rate + k_SO4 * mSO4/(1.e-4 + mSO4) +110 moles = rate * M * (M/M0) * TIME +200 SAVE moles + -end + +########### +#Pyrolusite +########### +# +# Postma, D. and Appelo, C.A.J., 2000, GCA, vol. 64, pp. 1237-1247. +# Rate equation given as mol L^-1 s^-1 +# +# Example of KINETICS data block for Pyrolusite +# KINETICS 1-12 +# Pyrolusite +# -tol 1.e-7 +# -m0 0.1 +# -m 0.1 +# -time 0.5 day in 10 +Pyrolusite + -start +10 if (M <= 0) THEN GOTO 200 +20 sr_pl = SR("Pyrolusite") +30 if (sr_pl > 1) THEN GOTO 100 +40 REM sr_pl <= 1, undersaturated +50 Fe_t = TOT("Fe(2)") +60 if Fe_t < 1e-8 then goto 200 +70 moles = 6.98e-5 * Fe_t * (M/M0)^0.67 * TIME * (1 - sr_pl) +80 GOTO 200 +100 REM sr_pl > 1, supersaturated +110 moles = 2e-3 * 6.98e-5 * (1 - sr_pl) * TIME +200 SAVE moles * SOLN_VOL + -end +END +# ============================================================================================= +#(a) means amorphous. (d) means disordered, or less crystalline. +#(14A) refers to 14 angstrom spacing of clay planes. FeS(ppt), +#precipitated, indicates an initial precipitate that is less crystalline. +#Zn(OH)2(e) indicates a specific crystal form, epsilon. +# ============================================================================================= +# For the reaction aA + bB = cC + dD, +# with delta_v = c*Vm(C) + d*Vm(D) - a*Vm(A) - b*Vm(B), +# PHREEQC adds the pressure term to log_k: -= delta_v * (P - 1) / (2.3RT). +# Vm(A) is volume of A, cm3/mol, P is pressure, atm, R is the gas constant, T is Kelvin. +# Gas-pressures and fugacity coefficients are calculated with Peng-Robinson's EOS. +# Binary interaction coefficients from Soreide and Whitson, 1992, FPE 77, 217 are +# hard-coded in calc_PR(): +# kij CH4 CO2 H2S N2 +# H2O 0.49 0.19 0.19 0.49 +# ============================================================================================= +# The molar volumes of solids are entered with +# -Vm vm cm3/mol +# vm is the molar volume, cm3/mol (default), but dm3/mol and m3/mol are permitted. +# Data for minerals' vm (= MW (g/mol) / rho (g/cm3)) are defined using rho from +# Deer, Howie and Zussman, The rock-forming minerals, Longman. +# -------------------- +# Temperature- and pressure-dependent volumina of aqueous species are calculated with a Redlich- +# type equation (cf. Redlich and Meyer, Chem. Rev. 64, 221), from parameters entered with +# -Vm a1 a2 a3 a4 W a0 i1 i2 i3 i4 +# The volume (cm3/mol) is +# Vm(T, pb, I) = 41.84 * (a1 * 0.1 + a2 * 100 / (2600 + pb) + a3 / (T - 228) + +# a4 * 1e4 / (2600 + pb) / (T - 228) - W * QBrn) +# + z^2 / 2 * Av * f(I^0.5) +# + (i1 + i2 / (T - 228) + i3 * (T - 228)) * I^i4 +# Volumina at I = 0 are obtained using supcrt92 formulas (Johnson et al., 1992, CG 18, 899). +# 41.84 transforms cal/bar/mol into cm3/mol. +# pb is pressure in bar. +# W * QBrn is the energy of solvation, calculated from W and the pressure dependence of the Born equation, +# W is fitted on measured solution densities. +# z is charge of the solute species. +# Av is the Debye-Hückel limiting slope (DH_AV in PHREEQC basic). +# a0 is the ion-size parameter in the extended Debye-Hückel equation: +# f(I^0.5) = I^0.5 / (1 + a0 * DH_B * I^0.5), +# a0 = -gamma x for cations, = 0 for anions. +# For details, consult ref. 1. +# +# ref. 1: Appelo, Parkhurst and Post, 2014. Geochim. Cosmochim. Acta 125, 49–67. +# ref. 2: Procedures from ref. 1 using data compiled by Laliberté, 2009, J. Chem. Eng. Data 54, 1725. +# ref. 3: Appelo, 2017, Cem. Concr. Res. 101, 102-113. +# +# ============================================================================================= +# It remains the responsibility of the user to check the calculated results, for example with +# measured solubilities as a function of (P, T). diff --git a/app/testEngine.pqi b/app/testEngine.pqi new file mode 100644 index 00000000..43e4eea7 --- /dev/null +++ b/app/testEngine.pqi @@ -0,0 +1,13 @@ +SOLUTION 1 +units mol/kgw +temp 25 +Ca 0.1 +Mg 0.1 +Cl 0.5 charge +Na 0.1 +PURE 1 +Calcite 0.0 1 +Dolomite 0.0 0 +## RUN_CELLS +## -cells 1 +END From 8c2e1ce8e86e0e9287f395cce5c58b286b8a76f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Sep 2024 15:18:49 +0200 Subject: [PATCH 4/6] refactor: Update printVectorsAsTable function in ipqEngine.cpp - Added "Array Index" column to the table output --- app/ipqEngine.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/app/ipqEngine.cpp b/app/ipqEngine.cpp index 7199ea92..e222b6d8 100644 --- a/app/ipqEngine.cpp +++ b/app/ipqEngine.cpp @@ -39,13 +39,14 @@ void printVectorsAsTable(const std::vector &names, return; } - std::cout << std::left << std::setw(20) << "Name" << std::setw(20) << "Value" - << std::endl; + std::cout << std::left << std::setw(20) << "Array Index" << std::setw(20) + << "Name" << std::setw(20) << "Value" << std::endl; std::cout << std::string(40, '-') << std::endl; for (std::size_t i = 0; i < names.size(); ++i) { - std::cout << std::left << std::setw(20) << names[i] << std::setw(20) - << values[i] << std::endl; + std::cout << std::left << std::setw(20) << std::to_string(i) + << std::setw(20) << names[i] << std::setw(20) << values[i] + << std::endl; } } @@ -111,8 +112,8 @@ int main(int argc, char *argv[]) { // std::for_each(NewInp.begin()+3, NewInp.end(), // [](double a) { return a * 0.8;}); - NewInp[7] += 0.2; // Cl - NewInp[9] += 0.2; // Na + NewInp[5] += 0.2; // Cl + NewInp[7] += 0.2; // Na MyEngine.runCell(NewInp, 200); std::cout << "\n\nAfter running the cell:\n"; From 9722151735e3a46acf1e0740bba46b0281a1fcad Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 10 Sep 2024 17:29:52 +0200 Subject: [PATCH 5/6] Selfcontained input in app/ipqEngine --- app/ipqEngine.cpp | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/app/ipqEngine.cpp b/app/ipqEngine.cpp index e222b6d8..e1bef11c 100644 --- a/app/ipqEngine.cpp +++ b/app/ipqEngine.cpp @@ -13,6 +13,7 @@ #include "PhreeqcEngine.hpp" #include "PhreeqcInit.hpp" + static std::string readFile(const std::string &path) { std::string string_rpath(PATH_MAX, '\0'); @@ -57,26 +58,42 @@ template void printVector(const std::vector &vec) { std::cout << std::endl; } + + +static std::string TestInput = R"(SOLUTION 1 + units mol/kgw + temp 25 + Ca 0.1 + Mg 0.1 + Cl 0.5 charge + Na 0.1 + PURE 1 + Calcite 0.0 1 + Dolomite 0.0 0 + ## RUN_CELLS + ## -cells 1 + END)"; + + int main(int argc, char *argv[]) { - if (argc < 3) { - std::cout << "Two args needed, PHREEQC script and database" << std::endl; + if (argc < 2) { + std::cout << "PHREEQC database needed as argument" << std::endl; return 1; } - const std::string &script = argv[1]; - const std::string &database = argv[2]; + const std::string &database = argv[1]; wordexp_t exp_result; - std::cout << "Reading script from file " << script << " using database " - << database << std::endl; + // std::cout << "Reading script from file " << script << " using database " + // << database << std::endl; // expand the script and database path - const std::string read_script = readFile(script); + // const std::string read_script = readFile(script); const std::string read_database = readFile(database); // initialize the module instance "ipqcmod" - PhreeqcInit ipqcmod(read_database, read_script); + PhreeqcInit ipqcmod(read_database, TestInput); // allocate tabular data PhreeqcInit::PhreeqcMat pqc_mat = ipqcmod.getPhreeqcMat(); @@ -87,7 +104,7 @@ int main(int argc, char *argv[]) { POETConfig MyConfig; MyConfig.database = read_database; - MyConfig.input_script = read_script; + MyConfig.input_script = TestInput; MyConfig.cell = POETInitCell{ // H, O and charge need to be removed. I'm not sure why I did it this way. std::vector(totals_names.begin() + 3, totals_names.end()), From e3c3e890df90d404096fb394db5a7b1856056ca1 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 10 Sep 2024 18:08:55 +0200 Subject: [PATCH 6/6] refactor: Enable std::map::contains only when cpp >= 20 is defined --- poet/include/PhreeqcInit.hpp | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/poet/include/PhreeqcInit.hpp b/poet/include/PhreeqcInit.hpp index 100992d4..02433678 100644 --- a/poet/include/PhreeqcInit.hpp +++ b/poet/include/PhreeqcInit.hpp @@ -102,7 +102,7 @@ public: * cell. Empty if no exchange is defined. */ std::vector getExchanger(int id) { - if (this->exchanger.contains(id)) { + if (contains_id(this->exchanger, id)) { return this->exchanger[id]; } @@ -117,7 +117,7 @@ public: * cell. Empty if no kinetics are defined. */ std::vector getKineticsNames(int id) { - if (this->kinetics.contains(id)) { + if (contains_id(this->kinetics, id)) { return this->kinetics[id]; } @@ -132,7 +132,7 @@ public: * cell. Empty if no equilibrium is defined. */ std::vector getEquilibriumNames(int id) { - if (this->equilibrium.contains(id)) { + if (contains_id(this->equilibrium, id)) { return this->equilibrium[id]; } @@ -147,7 +147,7 @@ public: * for the cell. Empty if no surface is defined. */ std::vector getSurfaceCompNames(int id) { - if (this->surface_comps.contains(id)) { + if (contains_id(this->surface_comps, id)) { return this->surface_comps[id]; } @@ -162,7 +162,7 @@ public: * the cell. Empty if no surface is defined. */ std::vector getSurfaceChargeNames(int id) { - if (this->surface_charge.contains(id)) { + if (contains_id(this->surface_charge, id)) { return this->surface_charge[id]; } @@ -216,6 +216,15 @@ private: std::map> surface_comps; std::map> surface_charge; + bool contains_id(const std::map> &map, + int id) const { +#if __cplusplus >= 202002L + return map.contains(id); +#else + return map.find(id) != map.end(); +#endif + } + std::set surface_primaries; using RawMap = std::map>>;