diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index a98eb8681fc5ebfb4c3f841d69758e2d38020772..448f87e594867f1e0c7d114e6cbab1c3bc6019be 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -74,8 +74,7 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
     // Find and extract mesh from the list of meshes.
     auto const& mesh = *BaseLib::findElementOrError(
         std::begin(meshes), std::end(meshes),
-        [&mesh_name](auto const& mesh)
-        {
+        [&mesh_name](auto const& mesh) {
             assert(mesh != nullptr);
             return mesh->getName() == mesh_name;
         },
@@ -110,6 +109,13 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
         output_directory,
         BaseLib::extractBaseNameWithoutExtension(config.getProjectFileName()));
 
+    if (!surface.empty() && !exchange.empty())
+    {
+        OGS_FATAL(
+            "Using surface and exchange reactions simultaneously is not "
+            "supported at the moment");
+    }
+
     auto dump = surface.empty() && exchange.empty()
                     ? nullptr
                     : std::make_unique<PhreeqcIOData::Dump>(project_file_name);
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 49e54b097245aa0cc6a61029e6bb8ca75b605f15..f48e3b393984f8a87baa23c73633313497fbe617 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -243,7 +243,7 @@ PhreeqcIO::PhreeqcIO(std::string const& project_file_name,
                      std::unique_ptr<ChemicalSystem>&& chemical_system,
                      std::vector<ReactionRate>&& reaction_rates,
                      std::vector<SurfaceSite>&& surface,
-                     std::vector<Exchange>&& exchange,
+                     std::vector<ExchangeSite>&& exchange,
                      std::unique_ptr<UserPunch>&& user_punch,
                      std::unique_ptr<Output>&& output,
                      std::unique_ptr<Dump>&& dump,
@@ -512,7 +512,7 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
         if (!exchange.empty())
         {
             os << "EXCHANGE " << chemical_system_id + 1 << "\n";
-            std::size_t aqueous_solution_id =
+            std::size_t const aqueous_solution_id =
                 dump->aqueous_solutions_prev.empty()
                     ? chemical_system_id + 1
                     : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
@@ -583,22 +583,15 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
     auto const& surface = phreeqc_io._surface;
     auto const& exchange = phreeqc_io._exchange;
 
-    int const reaction_case = !surface.empty() + !exchange.empty();
-    int num_skipped_lines;
-
-    switch (reaction_case)
+    if (!surface.empty() && !exchange.empty())
     {
-        case 0:
-            num_skipped_lines = 1;
-            break;
-        case 1:
-            num_skipped_lines = 2;
-            break;
-        case 2:
-            num_skipped_lines = 3;
-            break;
+        OGS_FATAL(
+            "Using surface and exchange reactions simultaneously is not "
+            "supported at the moment");
     }
 
+    int const num_skipped_lines = surface.empty() && exchange.empty() ? 1 : 2;
+
     auto& equilibrium_reactants =
         phreeqc_io._chemical_system->equilibrium_reactants;
     auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 7cfe9f48eadfbe920deb8445896f13254bca7df3..dc73f5fa1383bb044c936629d9dc589b80d8deb2 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -28,7 +28,7 @@ struct ChemicalSystem;
 struct ReactionRate;
 struct Output;
 struct SurfaceSite;
-struct Exchange;
+struct ExchangeSite;
 struct Dump;
 struct UserPunch;
 
@@ -40,7 +40,7 @@ public:
               std::unique_ptr<ChemicalSystem>&& chemical_system,
               std::vector<ReactionRate>&& reaction_rates,
               std::vector<SurfaceSite>&& surface,
-              std::vector<Exchange>&& exchange,
+              std::vector<ExchangeSite>&& exchange,
               std::unique_ptr<UserPunch>&& user_punch,
               std::unique_ptr<Output>&& output,
               std::unique_ptr<Dump>&& dump,
@@ -109,7 +109,7 @@ private:
     std::unique_ptr<ChemicalSystem> _chemical_system;
     std::vector<ReactionRate> const _reaction_rates;
     std::vector<SurfaceSite> const _surface;
-    std::vector<Exchange> const _exchange;
+    std::vector<ExchangeSite> const _exchange;
     std::unique_ptr<UserPunch> _user_punch;
     std::unique_ptr<Output> const _output;
     std::unique_ptr<Dump> const _dump;
diff --git a/ChemistryLib/PhreeqcIOData/CreateExchange.cpp b/ChemistryLib/PhreeqcIOData/CreateExchange.cpp
index 15bab28b168c3e48269062289c00c7362af93715..a51962314b3f151d051b2f1530c8218daf0e05b6 100644
--- a/ChemistryLib/PhreeqcIOData/CreateExchange.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateExchange.cpp
@@ -8,7 +8,7 @@
  *
  */
 
-#include <optional>
+#include "CreateExchange.h"
 
 #include "BaseLib/ConfigTree.h"
 #include "Exchange.h"
@@ -17,7 +17,7 @@ namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-std::vector<Exchange> createExchange(
+std::vector<ExchangeSite> createExchange(
     std::optional<BaseLib::ConfigTree> const& config)
 {
     if (!config)
@@ -25,20 +25,20 @@ std::vector<Exchange> createExchange(
         return {};
     }
 
-    std::vector<Exchange> exchange;
+    std::vector<ExchangeSite> exchange;
     for (auto const& site_config :
-         //! \ogs_file_param{prj__chemical_system__exchange_assemblage}
-         config->getConfigSubtreeList("exchange_assemblage"))
+         //! \ogs_file_param{prj__chemical_system__exchange__exchange_site}
+         config->getConfigSubtreeList("exchange_site"))
     {
-        //! \ogs_file_param{prj__chemical_system__exchange_assemblage__name}
-        auto name = site_config.getConfigParameter<std::string>("name");
+        //! \ogs_file_param{prj__chemical_system__exchange__exchange_site__ion_exchanging_species}
+        auto name = site_config.getConfigParameter<std::string>(
+            "ion_exchanging_species");
 
         auto const molality =
-            //! \ogs_file_param{prj__chemical_system__exchange_assemblage__molality}
+            //! \ogs_file_param{prj__chemical_system__exchange__exchange_site__molality}
             site_config.getConfigParameter<double>("molality");
 
-        exchange.emplace_back(
-            std::move(name), molality);
+        exchange.emplace_back(std::move(name), molality);
     }
 
     return exchange;
diff --git a/ChemistryLib/PhreeqcIOData/CreateExchange.h b/ChemistryLib/PhreeqcIOData/CreateExchange.h
index 447527cde1b452b447b03e211e1c6a2b5ebc6589..4fdbb9fc9a24342f6e564c9b6fd3790a974b3e59 100644
--- a/ChemistryLib/PhreeqcIOData/CreateExchange.h
+++ b/ChemistryLib/PhreeqcIOData/CreateExchange.h
@@ -22,9 +22,9 @@ namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-struct Exchange;
+struct ExchangeSite;
 
-std::vector<Exchange> createExchange(
+std::vector<ExchangeSite> createExchange(
     std::optional<BaseLib::ConfigTree> const& config);
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/Exchange.cpp b/ChemistryLib/PhreeqcIOData/Exchange.cpp
index 556132068919db7158c731d29bfb5ad64b651de3..ac7786785a31c0b4447805900ab5e232395d6b48 100644
--- a/ChemistryLib/PhreeqcIOData/Exchange.cpp
+++ b/ChemistryLib/PhreeqcIOData/Exchange.cpp
@@ -16,9 +16,9 @@ namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-std::ostream& operator<<(std::ostream& os, Exchange const& exchange_assemblage)
+std::ostream& operator<<(std::ostream& os, ExchangeSite const& exchange_site)
 {
-    os << exchange_assemblage.name << " " << exchange_assemblage.molality << "\n";
+    os << exchange_site.ion_exchanging_species << " " << exchange_site.molality << "\n";
 
     return os;
 }
diff --git a/ChemistryLib/PhreeqcIOData/Exchange.h b/ChemistryLib/PhreeqcIOData/Exchange.h
index f8eaa7008d7cdee992167475945111c1a7a1e500..ed591c95f0ef5890c64a1b5764e8d2be6d5035c7 100644
--- a/ChemistryLib/PhreeqcIOData/Exchange.h
+++ b/ChemistryLib/PhreeqcIOData/Exchange.h
@@ -22,18 +22,19 @@ namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-struct Exchange
+struct ExchangeSite
 {
-    Exchange(std::string name_, double const molality_)
-        : name(std::move(name_)), molality(molality_)
+    ExchangeSite(std::string ion_exchanging_species_, double const molality_)
+        : ion_exchanging_species(std::move(ion_exchanging_species_)),
+          molality(molality_)
     {
     }
 
-    friend std::ostream& operator<<(std::ostream& os,
-                                    Exchange const& exchange_assemblage);
-
-    std::string const name;
+    std::string const ion_exchanging_species;
     double const molality;
 };
+
+std::ostream& operator<<(std::ostream& os, ExchangeSite const& exchange_site);
+
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/i_exchange_site.md b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/i_exchange_site.md
new file mode 100644
index 0000000000000000000000000000000000000000..f9d1114bfe9601f9e3bd81b2b75d0ee5dd4dbca3
--- /dev/null
+++ b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/i_exchange_site.md
@@ -0,0 +1 @@
+Defines an exchange site.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_ion_exchanging_species.md b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_ion_exchanging_species.md
new file mode 100644
index 0000000000000000000000000000000000000000..bae20a715dfa5c30cfba9c61db1c36e1e80648fa
--- /dev/null
+++ b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_ion_exchanging_species.md
@@ -0,0 +1 @@
+This keyword is used to input the species of the exchange site. When using the Phreeqc chemical solver, the species should be defined in the EXCHANGE_MASTER_SPECIES keyword in the database.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_molality.md b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_molality.md
new file mode 100644
index 0000000000000000000000000000000000000000..5f1ab4186011bf93e625b9765d013763eb291399
--- /dev/null
+++ b/Documentation/ProjectFile/prj/chemical_system/exchange/exchange_site/t_molality.md
@@ -0,0 +1 @@
+This keyword is used to input the amount (mol/kg of water) of the exchange site.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/chemical_system/exchange/i_exchange.md b/Documentation/ProjectFile/prj/chemical_system/exchange/i_exchange.md
new file mode 100644
index 0000000000000000000000000000000000000000..f6dd4d82f5a22bd43dff7d5b242299bceb905830
--- /dev/null
+++ b/Documentation/ProjectFile/prj/chemical_system/exchange/i_exchange.md
@@ -0,0 +1 @@
+This keyword is used to define the amount and composition of a number of exchange sites.
\ No newline at end of file
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index 12c916635772f50c7faa9ccaa111fcdf158c8781..4322142a7ca3ee7e22c4bc7fe0650f74a8ad6489 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -1027,6 +1027,7 @@ if (NOT OGS_USE_MPI)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj RUNTIME 60)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/EquilibriumPhase/calciteDissolvePrecipitateOnly.prj RUNTIME 25)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/EquilibriumPhase/calcitePorosityChange.prj RUNTIME 25)
+    OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.prj RUNTIME 60)
 endif()
 
 AddTest(
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/ReactiveDomain_exchange.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/ReactiveDomain_exchange.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..01337c429947d28b2b064a4b7bedba8b0da826ad
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/ReactiveDomain_exchange.vtu
@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0"                    RangeMax="40"                   offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="448"                 />
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="0"                    RangeMax="39"                   offset="672"                 />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="1112"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="2436"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="3300"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="3740"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _SAEAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABMAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABYAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABkAAAAAAAAAGgAAAAAAAAAbAAAAAAAAABwAAAAAAAAAHQAAAAAAAAAeAAAAAAAAAB8AAAAAAAAAIAAAAAAAAAAhAAAAAAAAACIAAAAAAAAAIwAAAAAAAAAkAAAAAAAAACUAAAAAAAAAJgAAAAAAAAAnAAAAAAAAACgAAAAAAAAAoAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAQAEAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABMAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABYAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABkAAAAAAAAAGgAAAAAAAAAbAAAAAAAAABwAAAAAAAAAHQAAAAAAAAAeAAAAAAAAAB8AAAAAAAAAIAAAAAAAAAAhAAAAAAAAACIAAAAAAAAAIwAAAAAAAAAkAAAAAAAAACUAAAAAAAAAJgAAAAAAAAAnAAAAAAAAAA==2AMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB7FK5H4Xq0PwAAAAAAAAAAAAAAAAAAAAAuiPHSTWJgPwAAAAAAAAAAAAAAAAAAAAAqhvHSTWJwPwAAAAAAAAAAAAAAAAAAAACFS2q8dJN4PwAAAAAAAAAAAAAAAAAAAAB3ifHSTWKAPwAAAAAAAAAAAAAAAAAAAAAF7a1H4XqEPwAAAAAAAAAAAAAAAAAAAAB5SWq8dJOIPwAAAAAAAAAAAAAAAAAAAAB5pyYxCKyMPwAAAAAAAAAAAAAAAAAAAACngfHSTWKQPwAAAAAAAAAAAAAAAAAAAAAxsE+Nl26SPwAAAAAAAAAAAAAAAAAAAAAu361H4XqUPwAAAAAAAAAAAAAAAAAAAAArDgwCK4eWPwAAAAAAAAAAAAAAAAAAAACMPGq8dJOYPwAAAAAAAAAAAAAAAAAAAACgash2vp+aPwAAAAAAAAAAAAAAAAAAAADFmSYxCKycPwAAAAAAAAAAAAAAAAAAAACcyITrUbiePwAAAAAAAAAAAAAAAAAAAACme/HSTWKgPwAAAAAAAAAAAAAAAAAAAAAVkSCwcmihPwAAAAAAAAAAAAAAAAAAAAC4qU+Nl26iPwAAAAAAAAAAAAAAAAAAAABuw35qvHSjPwAAAAAAAAAAAAAAAAAAAADH2a1H4XqkPwAAAAAAAAAAAAAAAAAAAAA3+NwkBoGlPwAAAAAAAAAAAAAAAAAAAACFFAwCK4emPwAAAAAAAAAAAAAAAAAAAACEMDvfT42nPwAAAAAAAAAAAAAAAAAAAAC2T2q8dJOoPwAAAAAAAAAAAAAAAAAAAADpbpmZmZmpPwAAAAAAAAAAAAAAAAAAAACEi8h2vp+qPwAAAAAAAAAAAAAAAAAAAAA2p/dT46WrPwAAAAAAAAAAAAAAAAAAAABpxiYxCKysPwAAAAAAAAAAAAAAAAAAAACa5VUOLbKtPwAAAAAAAAAAAAAAAAAAAAD+AIXrUbiuPwAAAAAAAAAAAAAAAAAAAABhHLTIdr6vPwAAAAAAAAAAAAAAAAAAAADKnfHSTWKwPwAAAAAAAAAAAAAAAAAAAABkLYlBYOWwPwAAAAAAAAAAAAAAAAAAAAD9vCCwcmixPwAAAAAAAAAAAAAAAAAAAACXTLgeheuxPwAAAAAAAAAAAAAAAAAAAACY2U+Nl26yPwAAAAAAAAAAAAAAAAAAAACvZef7qfGyPwAAAAAAAAAAAAAAAAAAAABI9X5qvHSzPwAAAAAAAAAAAAAAAAAAAADhhBbZzvezPwAAAAAAAAAAAAAAAAAAAAA=gAIAAAAAAAAAAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAMAAAAAAAAABAAAAAAAAAAEAAAAAAAAAAUAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAYAAAAAAAAABwAAAAAAAAAHAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAkAAAAAAAAACgAAAAAAAAAKAAAAAAAAAAsAAAAAAAAACwAAAAAAAAAMAAAAAAAAAAwAAAAAAAAADQAAAAAAAAANAAAAAAAAAA4AAAAAAAAADgAAAAAAAAAPAAAAAAAAAA8AAAAAAAAAEAAAAAAAAAAQAAAAAAAAABEAAAAAAAAAEQAAAAAAAAASAAAAAAAAABIAAAAAAAAAEwAAAAAAAAATAAAAAAAAABQAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABUAAAAAAAAAFgAAAAAAAAAWAAAAAAAAABcAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABgAAAAAAAAAGQAAAAAAAAAZAAAAAAAAABoAAAAAAAAAGgAAAAAAAAAbAAAAAAAAABsAAAAAAAAAHAAAAAAAAAAcAAAAAAAAAB0AAAAAAAAAHQAAAAAAAAAeAAAAAAAAAB4AAAAAAAAAHwAAAAAAAAAfAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAhAAAAAAAAACEAAAAAAAAAIgAAAAAAAAAiAAAAAAAAACMAAAAAAAAAIwAAAAAAAAAkAAAAAAAAACQAAAAAAAAAJQAAAAAAAAAlAAAAAAAAACYAAAAAAAAAJgAAAAAAAAAnAAAAAAAAACcAAAAAAAAAKAAAAAAAAAAoAAAAAAAAAAEAAAAAAAAAQAEAAAAAAAACAAAAAAAAAAQAAAAAAAAABgAAAAAAAAAIAAAAAAAAAAoAAAAAAAAADAAAAAAAAAAOAAAAAAAAABAAAAAAAAAAEgAAAAAAAAAUAAAAAAAAABYAAAAAAAAAGAAAAAAAAAAaAAAAAAAAABwAAAAAAAAAHgAAAAAAAAAgAAAAAAAAACIAAAAAAAAAJAAAAAAAAAAmAAAAAAAAACgAAAAAAAAAKgAAAAAAAAAsAAAAAAAAAC4AAAAAAAAAMAAAAAAAAAAyAAAAAAAAADQAAAAAAAAANgAAAAAAAAA4AAAAAAAAADoAAAAAAAAAPAAAAAAAAAA+AAAAAAAAAEAAAAAAAAAAQgAAAAAAAABEAAAAAAAAAEYAAAAAAAAASAAAAAAAAABKAAAAAAAAAEwAAAAAAAAATgAAAAAAAABQAAAAAAAAAA==KAAAAAAAAAADAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMD
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.prj b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.prj
new file mode 100644
index 0000000000000000000000000000000000000000..eb2a96d0f7519357b0244fc111a08692519528cb
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.prj
@@ -0,0 +1,695 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>exchange.vtu</mesh>
+        <mesh>exchange_upstream.vtu</mesh>
+        <mesh>exchange_downstream.vtu</mesh>
+        <mesh>ReactiveDomain_exchange.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>hc</name>
+            <type>ComponentTransport</type>
+            <integration_order>2</integration_order>
+            <coupling_scheme>staggered</coupling_scheme>
+            <process_variables>
+                <concentration>Ca</concentration>
+                <concentration>Cl</concentration>
+                <concentration>H</concentration>
+                <concentration>K</concentration>
+                <concentration>N(5)</concentration>
+                <concentration>Na</concentration>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <specific_body_force>0</specific_body_force>
+            <secondary_variables>
+                <secondary_variable internal_name="darcy_velocity" output_name="darcy_velocity"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <media>
+        <medium id="0">
+            <phases>
+                <phase>
+                    <type>AqueousLiquid</type>
+                    <components>
+                        <component>
+                            <name>Ca</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>Cl</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>H</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>K</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>N(5)</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>Na</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Constant</type>
+                                    <value>0</value>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                    <properties>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value>1.00e+03</value>
+                        </property>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value>1.00e-03</value>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>permeability</name>
+                    <type>Constant</type>
+                    <value>1.73207e-12</value>
+                </property>
+                <property>
+                    <name>porosity</name>
+                    <type>Constant</type>
+                    <value>0.1</value>
+                </property>
+                <property>
+                    <name>longitudinal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0.002</value>
+                </property>
+                <property>
+                    <name>transversal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0.1</value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <time_loop>
+        <global_process_coupling>
+            <max_iter>6</max_iter>
+            <convergence_criteria>
+            <!-- convergence criterion for process p -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process Ca -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process Cl -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process H -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process K -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process N(5) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            <!-- convergence criterion for process Na -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+            </convergence_criteria>
+        </global_process_coupling>
+        <processes>
+            <!-- convergence criterion for process p -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process Ca -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process Cl -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process H -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process K -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process N(5) -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+            <!-- convergence criterion for process Na -->
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>72000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>400</repeat>
+                            <delta_t>180</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>exchangeOutput</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>4</repeat>
+                    <each_steps>100</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>Ca</variable>
+                <variable>Cl</variable>
+                <variable>H</variable>
+                <variable>K</variable>
+                <variable>N(5)</variable>
+                <variable>Na</variable>
+                <variable>pressure</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <chemical_system chemical_solver="Phreeqc">
+        <mesh>ReactiveDomain_exchange</mesh>
+        <database>phreeqc.dat</database>
+        <solution>
+            <temperature>25</temperature>
+            <pressure>1</pressure>
+            <pe>12.5</pe>
+            <components>
+                <component>Ca</component>
+                <component>Cl</component>
+                <component>K</component>
+                <component>N(5)</component>
+                <component>Na</component>
+            </components>
+            <charge_balance>pH</charge_balance>
+        </solution>
+		<exchange>
+			<exchange_site>
+				<ion_exchanging_species>X</ion_exchanging_species>
+				<molality>1.1e-3</molality>
+			</exchange_site>
+		</exchange>
+        <knobs>
+            <max_iter>500</max_iter>
+            <relative_convergence_tolerance>1e-12</relative_convergence_tolerance>
+            <tolerance>1e-15</tolerance>
+            <step_size>5</step_size>
+            <scaling>0</scaling>
+        </knobs>
+    </chemical_system>
+    <parameters>
+        <parameter>
+            <name>c0_Ca</name>
+            <type>Constant</type>
+            <value>1e-12</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_Ca</name>
+            <type>Constant</type>
+            <value>6.0e-04</value>
+        </parameter>
+        <parameter>
+            <name>c0_Cl</name>
+            <type>Constant</type>
+            <value>1e-12</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_Cl</name>
+            <type>Constant</type>
+            <value>1.20e-03</value>
+        </parameter>
+        <parameter>
+            <name>c0_H</name>
+            <type>Constant</type>
+            <value>1.0e-07</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_H</name>
+            <type>Constant</type>
+            <value>1.0e-07</value>
+        </parameter>
+        <parameter>
+            <name>c0_K</name>
+            <type>Constant</type>
+            <value>0.2e-3</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_K</name>
+            <type>Constant</type>
+            <value>1e-12</value>
+        </parameter>
+        <parameter>
+            <name>c0_N(5)</name>
+            <type>Constant</type>
+            <value>1.2e-3</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_N(5)</name>
+            <type>Constant</type>
+            <value>1e-12</value>
+        </parameter>
+        <parameter>
+            <name>c0_Na</name>
+            <type>Constant</type>
+            <value>1.0e-3</value>
+        </parameter>
+        <parameter>
+            <name>c_inlet_Na</name>
+            <type>Constant</type>
+            <value>1e-12</value>
+        </parameter>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>1e5</value>
+        </parameter>
+        <parameter>
+            <name>p_upstream</name>
+            <type>Constant</type>
+            <value>1e5</value>
+        </parameter>
+        <parameter>
+            <name>p_downstream_Neumann</name>
+            <type>Constant</type>
+            <!-- vel: porosity * cell length / time step from Phreeqc original benchmark-->
+            <value>-2.777777777777778e-4</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>Ca</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_Ca</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_Ca</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>Cl</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_Cl</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_Cl</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>H</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_H</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_H</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>K</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_K</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_K</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>N(5)</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_N(5)</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_N(5)</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>Na</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_Na</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_inlet_Na</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>exchange_upstream</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_upstream</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>exchange_downstream</mesh>
+                    <type>Neumann</type>
+                    <parameter>p_downstream_Neumann</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>100</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>100</max_iteration_step>
+                <error_tolerance>1e-14</error_tolerance>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>pressure</field>
+            <absolute_tolerance>1e-6</absolute_tolerance>
+            <relative_tolerance>1e-10</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>Ca</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>Cl</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>H</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>K</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>N(5)</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>exchangeOutput_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>Na</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
\ No newline at end of file
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..5d345424389766bedfc17409e129d7d69d69b688
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.vtu
@@ -0,0 +1,23 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="0"                   />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="224"                 />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1548"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="2412"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="2852"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _oAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA2AMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB7FK5H4Xq0PwAAAAAAAAAAAAAAAAAAAAAuiPHSTWJgPwAAAAAAAAAAAAAAAAAAAAAqhvHSTWJwPwAAAAAAAAAAAAAAAAAAAACFS2q8dJN4PwAAAAAAAAAAAAAAAAAAAAB3ifHSTWKAPwAAAAAAAAAAAAAAAAAAAAAF7a1H4XqEPwAAAAAAAAAAAAAAAAAAAAB5SWq8dJOIPwAAAAAAAAAAAAAAAAAAAAB5pyYxCKyMPwAAAAAAAAAAAAAAAAAAAACngfHSTWKQPwAAAAAAAAAAAAAAAAAAAAAxsE+Nl26SPwAAAAAAAAAAAAAAAAAAAAAu361H4XqUPwAAAAAAAAAAAAAAAAAAAAArDgwCK4eWPwAAAAAAAAAAAAAAAAAAAACMPGq8dJOYPwAAAAAAAAAAAAAAAAAAAACgash2vp+aPwAAAAAAAAAAAAAAAAAAAADFmSYxCKycPwAAAAAAAAAAAAAAAAAAAACcyITrUbiePwAAAAAAAAAAAAAAAAAAAACme/HSTWKgPwAAAAAAAAAAAAAAAAAAAAAVkSCwcmihPwAAAAAAAAAAAAAAAAAAAAC4qU+Nl26iPwAAAAAAAAAAAAAAAAAAAABuw35qvHSjPwAAAAAAAAAAAAAAAAAAAADH2a1H4XqkPwAAAAAAAAAAAAAAAAAAAAA3+NwkBoGlPwAAAAAAAAAAAAAAAAAAAACFFAwCK4emPwAAAAAAAAAAAAAAAAAAAACEMDvfT42nPwAAAAAAAAAAAAAAAAAAAAC2T2q8dJOoPwAAAAAAAAAAAAAAAAAAAADpbpmZmZmpPwAAAAAAAAAAAAAAAAAAAACEi8h2vp+qPwAAAAAAAAAAAAAAAAAAAAA2p/dT46WrPwAAAAAAAAAAAAAAAAAAAABpxiYxCKysPwAAAAAAAAAAAAAAAAAAAACa5VUOLbKtPwAAAAAAAAAAAAAAAAAAAAD+AIXrUbiuPwAAAAAAAAAAAAAAAAAAAABhHLTIdr6vPwAAAAAAAAAAAAAAAAAAAADKnfHSTWKwPwAAAAAAAAAAAAAAAAAAAABkLYlBYOWwPwAAAAAAAAAAAAAAAAAAAAD9vCCwcmixPwAAAAAAAAAAAAAAAAAAAACXTLgeheuxPwAAAAAAAAAAAAAAAAAAAACY2U+Nl26yPwAAAAAAAAAAAAAAAAAAAACvZef7qfGyPwAAAAAAAAAAAAAAAAAAAABI9X5qvHSzPwAAAAAAAAAAAAAAAAAAAADhhBbZzvezPwAAAAAAAAAAAAAAAAAAAAA=gAIAAAAAAAAAAAAAAAAAAAIAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAMAAAAAAAAABAAAAAAAAAAEAAAAAAAAAAUAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAYAAAAAAAAABwAAAAAAAAAHAAAAAAAAAAgAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAkAAAAAAAAACgAAAAAAAAAKAAAAAAAAAAsAAAAAAAAACwAAAAAAAAAMAAAAAAAAAAwAAAAAAAAADQAAAAAAAAANAAAAAAAAAA4AAAAAAAAADgAAAAAAAAAPAAAAAAAAAA8AAAAAAAAAEAAAAAAAAAAQAAAAAAAAABEAAAAAAAAAEQAAAAAAAAASAAAAAAAAABIAAAAAAAAAEwAAAAAAAAATAAAAAAAAABQAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABUAAAAAAAAAFgAAAAAAAAAWAAAAAAAAABcAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABgAAAAAAAAAGQAAAAAAAAAZAAAAAAAAABoAAAAAAAAAGgAAAAAAAAAbAAAAAAAAABsAAAAAAAAAHAAAAAAAAAAcAAAAAAAAAB0AAAAAAAAAHQAAAAAAAAAeAAAAAAAAAB4AAAAAAAAAHwAAAAAAAAAfAAAAAAAAACAAAAAAAAAAIAAAAAAAAAAhAAAAAAAAACEAAAAAAAAAIgAAAAAAAAAiAAAAAAAAACMAAAAAAAAAIwAAAAAAAAAkAAAAAAAAACQAAAAAAAAAJQAAAAAAAAAlAAAAAAAAACYAAAAAAAAAJgAAAAAAAAAnAAAAAAAAACcAAAAAAAAAKAAAAAAAAAAoAAAAAAAAAAEAAAAAAAAAQAEAAAAAAAACAAAAAAAAAAQAAAAAAAAABgAAAAAAAAAIAAAAAAAAAAoAAAAAAAAADAAAAAAAAAAOAAAAAAAAABAAAAAAAAAAEgAAAAAAAAAUAAAAAAAAABYAAAAAAAAAGAAAAAAAAAAaAAAAAAAAABwAAAAAAAAAHgAAAAAAAAAgAAAAAAAAACIAAAAAAAAAJAAAAAAAAAAmAAAAAAAAACgAAAAAAAAAKgAAAAAAAAAsAAAAAAAAAC4AAAAAAAAAMAAAAAAAAAAyAAAAAAAAADQAAAAAAAAANgAAAAAAAAA4AAAAAAAAADoAAAAAAAAAPAAAAAAAAAA+AAAAAAAAAEAAAAAAAAAAQgAAAAAAAABEAAAAAAAAAEYAAAAAAAAASAAAAAAAAABKAAAAAAAAAEwAAAAAAAAATgAAAAAAAABQAAAAAAAAAA==KAAAAAAAAAADAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMDAwMD
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_0_t_0.000000.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_0_t_0.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..d4e26fe2a6a47ecfb25a3ba2daa70c8c35e1a3c3
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_0_t_0.000000.vtu
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+      <DataArray type="Float64" Name="pe" NumberOfTuples="80" format="appended" RangeMin="12.559889828"         RangeMax="12.559889828"         offset="84"                  />
+    </FieldData>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="Float64" Name="Ca" format="appended" RangeMin="1.0000000487e-12"     RangeMax="1.0000000487e-12"     offset="160"                 />
+        <DataArray type="Float64" Name="Cl" format="appended" RangeMin="1e-12"                RangeMax="1e-12"                offset="232"                 />
+        <DataArray type="Float64" Name="H" format="appended" RangeMin="9.7789238893e-08"     RangeMax="9.7789238893e-08"     offset="304"                 />
+        <DataArray type="Float64" Name="K" format="appended" RangeMin="0.0002"               RangeMax="0.0002"               offset="376"                 />
+        <DataArray type="Float64" Name="N(5)" format="appended" RangeMin="0.0011999938369"      RangeMax="0.0011999938369"      offset="452"                 />
+        <DataArray type="Float64" Name="Na" format="appended" RangeMin="0.001"                RangeMax="0.001"                offset="524"                 />
+        <DataArray type="Float64" Name="darcy_velocity" format="appended" RangeMin="-6.8304736867e-18"    RangeMax="5.4701387734e-18"     offset="600"                 />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="100000"               RangeMax="100000"               offset="892"                 />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="964"                 />
+        <DataArray type="Float64" Name="velocity" format="appended" RangeMin="-6.8304736867e-18"    RangeMax="6.9253413768e-18"     offset="1024"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="1256"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1756"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="1920"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="2084"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9A1sTDQTTdPNUw2NUwzMzYHADNCBPA=AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAFwAAAAAAAAA=eF77MU3z4Uo5TYcfo/QoPQA0APxFLpA=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFQAAAAAAAAA=eF7jjO7tnzm90JZzlKaIBgCTRooPAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFQAAAAAAAAA=eF7L/KvbOHN6oW3mKE0RDQAt3qHDAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFQAAAAAAAAA=eF7rbL9WzeZQZdc5SlNEAwAzfYo4AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFwAAAAAAAAA=eF7LypN5/chMyz4TSmeN0mTRAB29imA=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFQAAAAAAAAA=eF5bHWQ1SXllsP3qUZoiGgA9fYFAAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFwAAAAAAAAA=eF5j2PPxkm9SgD0jlGYYpcmiAboLmUc=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAuQAAAAAAAAA=eF5jYACCGcF7GEE0w+fdYCokcg8TiPYw2wPmqwXYgGkGwz0cIMpBHCL+wwsiLqAAoS8YQuglgRA6JdQGbK6Kkw1YH8Nz63//gcCeHWKfiiNEnVUkhPYLgtALDPeA1YVbQfgeThDaKBRir1UYhD4gARGfYAPh10VDzBUJgvA1LCD0EzeIfyJs9iiABd5Yg9UJuNuAxV8YQNx5wh2i/ow/xBwgDebXQNVVGEHsYwuHiDfEw9UBAFa+O4k=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAEwAAAAAAAAA=eF5jYACCrB8ODKM0RTQAZ5xC8w==AQAAAAAAAAAAgAAAAAAAAKAAAAAAAAAADAAAAAAAAAA=eF5jYBjcAAAAoAABAQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAjQAAAAAAAAA=eF5jYACCGcF7QBSDQDyEfhNsA6YPeEDoByEQ8R+BED6DE4T+YAwRP+EN4aeEotIMWhD5BzoQvoMJhG8SBuXHQOgdDhBaJxAiHxEH4St4QPghcRDawQdCW0Dd0QB1b004VJ8nhGZgtkERXwB17wFHCP9DPKo8w+fdKObB+CfCIeokfPegyCOpAwAYVjQ4AQAAAAAAAAAAgAAAAAAAANgDAAAAAAAAVgEAAAAAAAA=eF5jYMAOqkXWuT+s2mKPLq7X8fGSb1IChrhWG0i8AEO81TtrT8nkCgzx8k6Q+gYMcda3a4H2tmCIV3qCzOnAFF+uZsixpgdDfHkjyPwJGOKGG/x7p+dNwvTXfZC9UzDEtfl4mLTbp2GI99iA3DMDQ3xB1omyffNnYYgfnQly5xwM8TknWl4H7piHIb6sGuT+BRjiohMVNhRlLMQQ37ES5K9FGOJ5h+uADl2MIX78Jsi/SzDEzX/cUWFrXIoZjyKgcFiGId5iYH3fv3c5hvg2f1D4rMAQf5k3EwhWYprTDQq3VRjiZsu/Bz9euhpDPPMYKDzXYIjPehrKp7tpLYb4P4ZWYDivwxBPlNkCtHg9hvipuaDw34AhnqLb6ZjwFFP87x5QvGzEEJ/us0Ou9TWm+IyboPjahCG+PvX575UfMcU9voLicTOG+MMWsZvnvmOKAwD0N598AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAWAAAAAAAAAA=eF5dxUkCgQAAAMAsoYUUZSkRLf7/QoecZi4TBLMlr3jNIW94yzuOOOaEU97zgTM+cs4Fn/jMJVd84Svf+M41N/zglp/84o7f/OGeBx554i8v/v8AsmsGaA==AQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAWgAAAAAAAAA=eF4txdEGwgAAAMAkk8lMkklmZmaSZGaSyczMJNn/f00P3b3cevW3ceCtQ+8cOfbeBx+d+OSzU2fOXbh05Yuvvvnu2o1bP/x055d7Dx49efbbH3+9+AfO3gZpAQAAAAAAAAAAgAAAAAAAACgAAAAAAAAADAAAAAAAAAA=eF5jZiYOAAAJxAB5
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_100_t_18000.000000.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_100_t_18000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f08895ae295a5931ec9b450e846452582a2406f6
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_100_t_18000.000000.vtu
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+      <DataArray type="Float64" Name="pe" NumberOfTuples="80" format="appended" RangeMin="12.562472773"         RangeMax="12.645428481"         offset="84"                  />
+    </FieldData>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="Float64" Name="Ca" format="appended" RangeMin="-6.9503331608e-07"    RangeMax="0.00060000000001"     offset="888"                 />
+        <DataArray type="Float64" Name="Cl" format="appended" RangeMin="4.5984181669e-05"     RangeMax="0.0012"               offset="1384"                />
+        <DataArray type="Float64" Name="H" format="appended" RangeMin="8.373178294e-08"      RangeMax="1.0117767633e-07"     offset="1880"                />
+        <DataArray type="Float64" Name="K" format="appended" RangeMin="9.7690317565e-13"     RangeMax="0.0007844158993"      offset="2376"                />
+        <DataArray type="Float64" Name="N(5)" format="appended" RangeMin="1.0000004589e-12"     RangeMax="0.0011539770009"      offset="2872"                />
+        <DataArray type="Float64" Name="Na" format="appended" RangeMin="9.9658327776e-13"     RangeMax="0.00099999975675"     offset="3368"                />
+        <DataArray type="Float64" Name="darcy_velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="3860"                />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="99987.170136"         RangeMax="100000"               offset="4152"                />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="4588"                />
+        <DataArray type="Float64" Name="velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="4648"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="4936"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="5436"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="5600"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="5764"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9A1sTDQTTdPNUw2NUwzMzYHADNCBPA=AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAOAIAAAAAAAA=eF4V0V1I02EYBXAIw82ZaBil4abzol7bYBpprNTXrzVDmyJZpDMslzNcMMuUXItQwxYRI5w1usmJbSupxoRMKkbOD9pyU1GSSWaDQJZTcZU0s/O/+nHO1eF5uiy1Z8cFhPIT66+EeYSye30bm1mE1rdf3nJmEBoryZXIswl9LZ85/0xMqIebz9FTQi8oxW/8OYReG2uZs+UT2h3htj7II/SnurnhaCGhDr1wVVFA6G1t5eeqIkI7Oby4d+g5AlN0toTQFb3c4kC/ZfjXGUR+Elg7fR+KbY7Ht04Qqj7C3SeCy7yY/hCcmTpQY4dX28penpESmiaKrE6B4d1xPCv86kmo0EDlzuGuNegSV7VOwOBDAet4MaHSSrcgBg5pDP3bMM/SqDPADN6oqfwkoe2a0prD8Cm/INgIX93Y3lML73Hjjc3wg6bgUhMUrWv9ahi62dPA9P42T89FeMhZuEsFZfHahBKodEnM5+DYyEKpEBZr7TO5kJ+0Q8WGG6IoIxcGgrLpJexw6vpKNuGYzjwwDL8nejOnINfkme9mdtvGG5/DNFXUUBNUBD49ugNjpB3JFXDwd2CxDo78aWnPgqkdZE0CjQP7V3lwPbmaLYR15nRxNCz6O5q6Fzpibawwc0dzUigCOZb1PrSC3FL6LeUXc1en664fygbf5izDHwqv3Ad96SrPIjzV+SJylvnHx97WL/C6Y/KYl/lDWfnUNLTXLk24mTyXaZyE2VbdQRd0z88uMP4HCjsDXA==AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/oMPVzAqqUM/b4T8+pqXcT0AwJ5qQKhDP3IzGhYrpkM/DBbXc+ahQz/+jF0XvJlDP7qqtnfcikM/5uTZgMVwQz9U7bcGXkRDP1P2OXej+kI/6gIJy8CCQj8EcC6dQMNBP23RPyM1lUA/zOJc1ZJ8PT+jMm434fU3P8bd3ogIvC8/fco+s/Cq/z6lIepsSlKnvjDrCAlxoGM+SKFbpCuYNz67Dxyi05kKPidAjlG4cN893s2x0cK1tD2EXxICJ9CRPXzUPBmc9ns9UFQQXrkSdD0waeHsB2RyPUNYjYah7HE9Ns5O4SPAcT2CSswA6atxPS/X8q7hoXE9cvqfgMiccT36Mb9lMJpxPS1fAx3hmHE9h1sFpjmYcT2nCfD55pdxPQj+yKW+l3E9Qn09NauXcT1iYLfyoZdxPTRBY5edl3E9qtHcl5uXcT3JQ5moAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/oE8VTAqqVM/K6+Dc+QbCD8ys8KaBKlTP2+dqCuwqFM/BMqdVQeoUz+yuP+0zqZTP/bB01ispFM/0FhqcB2hUz9aeZbsaptTP+4CEhaeklM/mvUllHaFUz92JjHJY3JTP9VzF6uDV1M/3lQAOKkyUz+SknVFbAFTP5blw55DwVI/uMRwP6lvUj8SbfUMRwpSP9sGd/goj1E/P78bJPD8UD+rbAzIAVNQP4oqT65YI08/G+bT0nt0TT8I8CJyEJ5LPyaaysegpkk/TG0JxVWWRz/0nsrfjHZFP5QhpydVUUM/LkNmwt8wQT9mtkIu4j0+PwZKbMK+SDo/HGvfmVKRNj8oSSQHWCMzP8aU9m74BjA/EqTdy4eBKj/P1wCEz6MlPwSFkF43cSE/IgBN5x/DGz/GKhtQjNQVP5RfqCr6BxE/rHf13Df3Cj/fUoqxAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/mZNZYHgKHs+rzxyGgF6dj5D/DA52Ch7PnJlR3XFKHs+ekmRc58oez4CCdrcVyh7Plia2hPYJ3s+tGtZh/0mez50P0T6lCV7PvYSkIRUI3s+sKMAUNQfez4viVrihRp7PkBqjsCoEns+BlrQZD4Hez60OwaQGfd6Prz04+7G4Ho+rJezro7Aej5CLVWofaR6PmcbwKzUg3o+bD65D05cej6SqF8+0y16PurXvHKm+Hk+OO92inS9eT7BimUgTX15PtI1TPmNOXk+4a38U8XzeD7k2GWJkK14PiZSknd8aHg+DgL00+sleD5cmROBBed3PnQXUjerrHc+pFLte3d3dz5A0AsxwUd3PrA7PPSiHXc+Akt8vAT5dj6FE3hpptl2PiivQYcqv3Y+3lrkfCGpdj5Y4RoJGJd2PlqoWb+3iHY+YLhJozV+dj6QspNqAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/oyL262UL3E974VBAuU2Kj9uCeP7xRiBPk5Sw8goIJw+DquWyUkesT74M5WdwkjCPoIy+ILyFdI+Tjys10v64D46E+TMv6nuPgC83D0r3fo+onVU1iT2Bj9yuNs/NDkTP9DcJcCJqh8/KhkgUSe6KT9WoZMWuV80P+ygaGfmZT8/fT0lUSi0ST+a6aH+efBHP8DB62EkLkQ/wg/MuFOOQD/WdcH6kLk6PyZDD2WIvjU/dUsgIFZJMj+UB0HQjxcwP7UscFyljy0/y9/m9ocNLD+l9pRgODUrP/Ip6IVevio/GXQp0hp+Kj+putFh2VsqP25trDrVSSo/UqkXzXpAKj+WDvptrzsqP+4FcyFCOSo/X5wkgQs4Kj96lHEYcjcqP4LMPj4nNyo/fj18JwM3Kj8o1Fj18TYqPz5mI93pNio/dcsXJuY2Kj+t74yYAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/lvsoQial3E94ZOoXiHoUj88+w3Oz8piPqrATwcOgX4+jJgN05Mtkj4WdUt+v9uiPnYZMclC97E+LovbnGYZwD6CsflDW37LPnl/7hn1i9Y+vHJs86/Z4T4IewkqBGPrPvK1HQ99afQ+tp76Dwig/T50uxkWlvcEP4c8nmya/Aw/tPy0WeiXEz+/n3Wr/e0ZPz6349jmzyA/m7WP/6JhJT+Ktc12CbEqP55AflfTXTA/9nqVBYW7Mz/Otrb7Umg3PyT9Z+0oVzs//NE1K7V3Pz+a1F6AnttBP2gDR0PRAEQ/2Ou30kEhRj8qoZfeKzNIP17vWz65LUo/lE7bU2sJTD/6QCb/ZMBNP7KwopKRTk8/b89ZgtRYUD90G9BTivRQP4Ni9k/celE/sjEoYNDsUT+74tXjuEtSP2S20XSBmFI/+NYO7kbRUj9tPZ8tAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUAEAAAAAAAA=eF7jDzriZdZRaBtlJzLBNynAPuBt0lUWhwq7dd4PNzdfnGwXt2rlkZCHy+34D5lff6C2084g9GtAb/4Ju3kCFyMkCq7Z/VPV8gjd89juz3cVneYjH+xSxSad5F/3226L39ZvJ73Y7Hum/pA8PkXQXmzhxwN7LaXt17lr1URtUrFXTK9YVsyqb39vzzmuX27m9jPVVgmw37azv3DE588kZWf7Q/UnBV4edrM/N0ltYsdXT/tVTBUb+Yt97IO2/7VZqetnP83BoLTHzd++o/IVR/Y7f/tYmTTDSvUA+5ZEO5VmpwB7H68Z8XGBAfYZIkLX0iID7CvM5FjWxwbY/9zA3Pk+PsB+8ctt16QSA+z5FGc7bQHSynk6bxiB/j823eyABpB+6+H71hpIVz2d7OEKpO/l5Zz2BNI/OTd+9AbRd7Z/8QHSXZkzs0HhBgAYBI+YAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAugAAAAAAAAA=eF4TDdWZpr9kkp1wjByYdvqiBaaX7dIF0x+MjMC0ngCEfvldH0z/y4Dwm1eYQfTNMAfTDKctwHRKtSWYDomD0Nt6IPSeVdZgui3bCkyvkYGoz86Eym+G0GbREHMffDEB012iEHUnTkD076qG2LcuxRBM26TpgekF7hDaoAri/m2MEP8ZXYK4+44NhE6ygMiHnIDQscsgdPRqiP9dfdXAdBqUXrcLQh+KUQHT804qgellcgpgGgDp9YuIAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAJAEAAAAAAAA=eF5jYACCrB8OrNsf7DDK/OGQeNHi3i8gfXhR4Z6vQFq5ZNWsD0C62fVxxWsg/VhMJvwZkHZ5EWL6EEgv3tkrfAdIM3Uf+3gNSCfE/D9/EUjv07VYdwZIy/wv6DkOpKsurMw6BKRvLHzksRdImxZLq+8A0pNdQlg3A+mPor2P1wJp/+dHD64A0mt3/Ju/GEjzdJnXzQPS2dEFMTOB9EmdlVZTgLTav4cS/UC67bzU904g/WxB8NUWIO1e1LO5Hkgvdz46sQpIc4r+yy8F0lnPzPwKgPSp7fk62UBau3MFVxqQ7ol6+CIBSL/VljoeDaT9/gYtDQPSG851NwcCaeEFRxJ9gHRp4V97d5D7nczknIC0jUj+HxsgveDp8lvmQBoAcQ2tMQ==AQAAAAAAAAAAgAAAAAAAAKAAAAAAAAAADAAAAAAAAAA=eF5jYBjcAAAAoAABAQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAtwAAAAAAAAA=eF4TDdWZpr9kkl3+ZE0wffuhIZhmbjEC0xFzIPyqYD0w3VxnBqEvQWjGOAsw/c/cEkyf2Aeh7/2FiPdLW4Hpl5a2YFp5rjmYNpkDkV9kBpG/bgDRN7XVGExrp0DMf34MIr7iFET/GR2IeyZA3fVMXxvizngIn2MqxB/JeRD3MplBxLc764JpXl0IPTsZIi//EuL/6bkaYPqTjgqYvg3lN7KrgmmWTcpg2vKhIpgWipED0wBu+IY5AQAAAAAAAAAAgAAAAAAAANgDAAAAAAAAVgEAAAAAAAA=eF5jYMAOqkXWuT+s2mKPLq7X8fGSb1IChrhWG0i8AEO81TtrT8nkCgzx8k6Q+gYMcda3a4H2tmCIV3qCzOnAFF+uZsixpgdDfHkjyPwJGOKGG/x7p+dNwvTXfZC9UzDEtfl4mLTbp2GI99iA3DMDQ3xB1omyffNnYYgfnQly5xwM8TknWl4H7piHIb6sGuT+BRjiohMVNhRlLMQQ37ES5K9FGOJ5h+uADl2MIX78Jsi/SzDEzX/cUWFrXIoZjyKgcFiGId5iYH3fv3c5hvg2f1D4rMAQf5k3EwhWYprTDQq3VRjiZsu/Bz9euhpDPPMYKDzXYIjPehrKp7tpLYb4P4ZWYDivwxBPlNkCtHg9hvipuaDw34AhnqLb6ZjwFFP87x5QvGzEEJ/us0Ou9TWm+IyboPjahCG+PvX575UfMcU9voLicTOG+MMWsZvnvmOKAwD0N598AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAWAAAAAAAAAA=eF5dxUkCgQAAAMAsoYUUZSkRLf7/QoecZi4TBLMlr3jNIW94yzuOOOaEU97zgTM+cs4Fn/jMJVd84Svf+M41N/zglp/84o7f/OGeBx554i8v/v8AsmsGaA==AQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAWgAAAAAAAAA=eF4txdEGwgAAAMAkk8lMkklmZmaSZGaSyczMJNn/f00P3b3cevW3ceCtQ+8cOfbeBx+d+OSzU2fOXbh05Yuvvvnu2o1bP/x055d7Dx49efbbH3+9+AfO3gZpAQAAAAAAAAAAgAAAAAAAACgAAAAAAAAADAAAAAAAAAA=eF5jZiYOAAAJxAB5
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_200_t_36000.000000.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_200_t_36000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f3bdb5f3d3222ea80c94ee494b04c58733192e4d
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_200_t_36000.000000.vtu
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+      <DataArray type="Float64" Name="pe" NumberOfTuples="80" format="appended" RangeMin="12.562472098"         RangeMax="12.646689571"         offset="84"                  />
+    </FieldData>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="Float64" Name="Ca" format="appended" RangeMin="-2.3660381338e-06"    RangeMax="0.0006"               offset="884"                 />
+        <DataArray type="Float64" Name="Cl" format="appended" RangeMin="0.0010539551624"      RangeMax="0.0012"               offset="1376"                />
+        <DataArray type="Float64" Name="H" format="appended" RangeMin="9.7664939221e-08"     RangeMax="1.0117767633e-07"     offset="1848"                />
+        <DataArray type="Float64" Name="K" format="appended" RangeMin="9.9999993374e-13"     RangeMax="0.000980479658"       offset="2312"                />
+        <DataArray type="Float64" Name="N(5)" format="appended" RangeMin="9.9999999995e-13"     RangeMax="0.00014603840963"     offset="2808"                />
+        <DataArray type="Float64" Name="Na" format="appended" RangeMin="9.9999999472e-13"     RangeMax="0.00095867964023"     offset="3304"                />
+        <DataArray type="Float64" Name="darcy_velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="3800"                />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="99987.170136"         RangeMax="100000"               offset="4092"                />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="4528"                />
+        <DataArray type="Float64" Name="velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="4588"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="4876"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="5376"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="5540"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="5704"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9A1sTDQTTdPNUw2NUwzMzYHADNCBPA=AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAANwIAAAAAAAA=eF4VzWtIE1AYxvGlTSkbtjScS1tsSR2ZmeVt2fLMLqaCgpVTzPDSRlhOSoicNCoWkq3StrFim7DSZDYiBkrWalguJbKB+SFNMkRdItvEUbB1sWeffvwfOO+5JR8sGhMSOpQ6Kv3DI7TUy/IGcwitGDb6XfsItSXbxDViQu1Gr7TvAKHThqaILkooCRmfzx8ilPWU6bQXEHolxme8JyH0jXw1L/cIoSlpuknZYULZTYL26qOEdp/8OevAXnfjouXgMUKjO1SVw9h7qr5/9KG79OnjGlg9e+6qqhD3euZE6VBNO5cDsFI19tsOddmNjFPHCTVIsgt4sMad+rUXZoobtrbCHzFW8xLM25bFHoEW6X0bvwjvees4kTBf2eMuh2XcZmEu1L7b1d4GXRwdUwY78x95zXCG82mhAxbKHP6XkGpbPf3w15asus9w59q/BRe0xjECHnjeazHNwEAJlxuEZy1Wmx9WKZaU0cWEhh6GDGvoOy+yjHHovS21K5ugf3wtbTssdgoqOJDU6yJ2w3zBVDMfvpILG/bD01Ga60L4IeObqAjOnTDMiaCZPzrkholMdtTlsPp47eYSQk2W1fq/aMWFUisX7VQnKDlwRbOhMBFmL4oWwvt0eXJbfHhPyrmbAK8pYs0sODD1bIINa59E2NfDdN+lpI3wiyCyMYi7i2ccmQx0ys1J2TI61v5+IhD+r2XENA133O5/MA+7yzJ6x6BMrvZPQI9kUP8a9u0ZePwW/gd8K/uUAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUAEAAAAAAAA=eF7rMQo10FrpbL/MxmdN1vIS22OX0i9qAvlSkie/aADpT5++KasD6YqFOV3KQHpG2rTNMkB6c+YGOQEg/cTnZv/PFc72RmIcPy4Caemz53b0Aul9PIYvJYH0ZJvtAhHLne2nHD3ZIr7M2d5vfaQq9xJn+1l5Kb07Fzjbf08vPpMy29neQbnb+ukkZ3ueWXI7C1qd7Y2PfdF1KXC235VkJ6Hg72x/qEJWv1Ha2b5tRcq0ktNO9hZZS7/e9neyV30/qfzkREf7W18qdgbkOdj/8/yVFpNqa5/BlPNvsa2FvcmM9gb5WgN79qmpeZ+c+Ow39fpbsN88vO8Um9HMy9Mq7Qyi3P+sFAiy67jv/e1DpZqd5ISj+SGff9mWBOywzrS7YCuxI9mB59Yy244jHedW/euyXcRnzix/sM5Wxecnl51nhS2btuWeS12ltgDbvZATAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAQQEAAAAAAAA=eF5TMAw10FoZbL/pE8/e2S6B9rlP12mA+M/+cYuD6AzvyE+aQPpi2bFVIDqBtV0VRGtO8DTTANIR3xKmqgHpv5z2F5WB9IclOZbyQFrGeM0BcSCt0hF5mRtIx5+slfu5Itg+5GGe8x0gHRTmXrwBSP/fpfQ2D0iXCH6O5APStwJFGJuWA/XbLLm4Y1mwfXK4vd2cpcH2e29pvuZbEmy/652r0euFwfZv0ypaGOcH2ycJWG4Vmx1sz5173unR1GD7U7+vbObpD7aPm7VIQ68t2N426ILL2apg+x6NgzPjcoLt5VZf3O8SFWy/u9Yp+KIzUF2Fl5y7RrD9xMzz3++yBdsHSH+cdvJekP20rRGr3NYH2dt5lxluLw+yFzWJzn5mHmR/JXy789YPgfYKexWP3lwaaO/jdWRXVlKgPQATSI35AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAOgEAAAAAAAA=eF67Mzm18YFGtZ1I08Pj1eZVdmZWLvUgvtZC60oQ7SWTkw+i13sERIPoxfeNTEH0a71tH+4D6Wt3KstBNFPmx6X3gHSLr43rXSCd8+jFz1tAWoZD3+Y6kJ4kP83iIkg8w3bdcSA920h47w4g/X1G4uyFQLpJykGnDmTv3qYcDyB9+GfqnO/q1XYci5vZ6oF0Xe2niDNq1XYr7iqbHFWttksyk9jurFJtNysiPltQqdpux+bY/PtyQPkdGf12UtV2JqpTpiwQBuqT8JktwlVt1xcmNDPsd5XdwiTWDa2fquzOsN++vuh1lZ3I0WMlZx9V2U3fkaZz5XqV3cojN8RiTlXZPXn7v59nV5Wd3KNFmrOXVdnZbtk98Xl/lV22ZK+YfVmVncH39QdYo6vsPonptyc5VdkBAOAOkic=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/vkpn22Zl3E93Tp2Gl6hLz/KZeNQgy7xPdbAsbBjgAw+6mQirfihIT7AsbX+DUIzPo5hbdcNlEM+PCNxukH5Uj78u8t9KMFhPqwFLCG4KnA+FteSTTfMfD5MIO9zbC2JPnarq0d7qpU+/rLBoY9joj7WueWrxNauPjKO4/Mrlbk+6HokWSUFxT6gbvT9uR/RPo0HVc62r9s+tACy7Fs65j5wZ1XgIrzxPkbNAFweJPw+bk6npy42Bj9r+YgJMnIRP6ywc2WRRxs/UCZ2NRk+JT+yP9JFWnMwP0Q+i/FOKDk/lii1liYpQz96JO5u+CNOP06NMQRuEFA/pOZL6ZbaTT+2PoE7q2NLP4YYjttSl0g/wc2fNweMRT8b/D6+OWpCP6yTbAdG1j4/c0bYC0KYOT9266uRhXA1P35U+guScjI/N57nQDyGMD8ac5KVAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/h/2KYGZl3E9Aa/9YjwkIz8SJ7XbbrK/PVzmfNZXj9k9HOfA4u8g7z0eOzUChroAPo6LwpkrsRA+MXLWr+qkHz5qnwVXo9ksPkwvhE1GfTk+KLHf2TfuRT7Zgssvz29SPgxu4JEMXV4+mokyPuOIaD7Sz24iKHtzPubydbQdb34+pyPaAkRohz6irKhwSb2RPgD1YVHGgpo+vsbXOdKKoz6A04p2W2+sPlIOzFDtbLQ+kPbY7777vD5UAImZLlHEPo5lJfG7Jsw+FG8cucRH0z78QFUFRB3aPkLLcC7dfeE+QtvUX0Eu5z7G+toNhWXuPi6rP4YtufM+zqdMCvpV+T69pnJvOhwAPw4ntclTSQQ/zmvBztFLCT8uXJZIHjwPP58opHr8FhM/qNmI5PwTFz9wiFHLAYcbPyB8EJ9VGiA/3fa5X8c1Ij8eb5acAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/oPOnn+Zl3E9L+be2fxpTz+Wzl/rEMDaPdInToBuBfY9G9m03HwVCz6aoVAuHV8dPq5cEO8kmi0+4LCqQb9jPD6aa8Zcbz1KPnQ3GF/5jVc+VuX7/ZqjZD6rM/NWE7dxPrSNupYd3n0+gRC6hETHiD5OIRwxL0OUPkBMvQfEWqA+1ENMFI4Wqj4u+1tm7ZS0PqFYveuBE8A+o+sSKZbiyD5I++8ijRnTPqKeHFWIGN0+mP5QoBIC5j7GRBypxorwPkF8CL7oufg+AKhNF8piAj9Wlgd+rTQLP/JljVRHBBQ/vIUtVHJgHT9QIBmU9pclP4L3OJ+pZC0/u7R2FkjsMj8yky/Hwdw3PzxNdZLsdT0//pbsekvGQT9yzoliGuhEPx4DoFIx50c/PkB4WDOGSj+KCQOXEZpMP0QsRFoLGU4/hCD6PzYPTz/vao9TAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAugAAAAAAAAA=eF4TDdWZpr9kkp1wjByYdvqiBaaX7dIF0x+MjMC0ngCEfvldH0z/y4Dwm1eYQfTNMAfTDKctwHRKtSWYDomD0Nt6IPSeVdZgui3bCkyvkYGoz86Eym+G0GbREHMffDEB012iEHUnTkD076qG2LcuxRBM26TpgekF7hDaoAri/m2MEP8ZXYK4+44NhE6ygMiHnIDQscsgdPRqiP9dfdXAdBqUXrcLQh+KUQHT804qgellcgpgGgDp9YuIAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAJAEAAAAAAAA=eF5jYACCrB8OrNsf7DDK/OGQeNHi3i8gfXhR4Z6vQFq5ZNWsD0C62fVxxWsg/VhMJvwZkHZ5EWL6EEgv3tkrfAdIM3Uf+3gNSCfE/D9/EUjv07VYdwZIy/wv6DkOpKsurMw6BKRvLHzksRdImxZLq+8A0pNdQlg3A+mPor2P1wJp/+dHD64A0mt3/Ju/GEjzdJnXzQPS2dEFMTOB9EmdlVZTgLTav4cS/UC67bzU904g/WxB8NUWIO1e1LO5Hkgvdz46sQpIc4r+yy8F0lnPzPwKgPSp7fk62UBau3MFVxqQ7ol6+CIBSL/VljoeDaT9/gYtDQPSG851NwcCaeEFRxJ9gHRp4V97d5D7nczknIC0jUj+HxsgveDp8lvmQBoAcQ2tMQ==AQAAAAAAAAAAgAAAAAAAAKAAAAAAAAAADAAAAAAAAAA=eF5jYBjcAAAAoAABAQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAtwAAAAAAAAA=eF4TDdWZpr9kkl3+ZE0wffuhIZhmbjEC0xFzIPyqYD0w3VxnBqEvQWjGOAsw/c/cEkyf2Aeh7/2FiPdLW4Hpl5a2YFp5rjmYNpkDkV9kBpG/bgDRN7XVGExrp0DMf34MIr7iFET/GR2IeyZA3fVMXxvizngIn2MqxB/JeRD3MplBxLc764JpXl0IPTsZIi//EuL/6bkaYPqTjgqYvg3lN7KrgmmWTcpg2vKhIpgWipED0wBu+IY5AQAAAAAAAAAAgAAAAAAAANgDAAAAAAAAVgEAAAAAAAA=eF5jYMAOqkXWuT+s2mKPLq7X8fGSb1IChrhWG0i8AEO81TtrT8nkCgzx8k6Q+gYMcda3a4H2tmCIV3qCzOnAFF+uZsixpgdDfHkjyPwJGOKGG/x7p+dNwvTXfZC9UzDEtfl4mLTbp2GI99iA3DMDQ3xB1omyffNnYYgfnQly5xwM8TknWl4H7piHIb6sGuT+BRjiohMVNhRlLMQQ37ES5K9FGOJ5h+uADl2MIX78Jsi/SzDEzX/cUWFrXIoZjyKgcFiGId5iYH3fv3c5hvg2f1D4rMAQf5k3EwhWYprTDQq3VRjiZsu/Bz9euhpDPPMYKDzXYIjPehrKp7tpLYb4P4ZWYDivwxBPlNkCtHg9hvipuaDw34AhnqLb6ZjwFFP87x5QvGzEEJ/us0Ou9TWm+IyboPjahCG+PvX575UfMcU9voLicTOG+MMWsZvnvmOKAwD0N598AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAWAAAAAAAAAA=eF5dxUkCgQAAAMAsoYUUZSkRLf7/QoecZi4TBLMlr3jNIW94yzuOOOaEU97zgTM+cs4Fn/jMJVd84Svf+M41N/zglp/84o7f/OGeBx554i8v/v8AsmsGaA==AQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAWgAAAAAAAAA=eF4txdEGwgAAAMAkk8lMkklmZmaSZGaSyczMJNn/f00P3b3cevW3ceCtQ+8cOfbeBx+d+OSzU2fOXbh05Yuvvvnu2o1bP/x055d7Dx49efbbH3+9+AfO3gZpAQAAAAAAAAAAgAAAAAAAACgAAAAAAAAADAAAAAAAAAA=eF5jZiYOAAAJxAB5
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_300_t_54000.000000.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_300_t_54000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..be9e3da8dbd439da054ed257b459e172b9ef6fe4
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_300_t_54000.000000.vtu
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+      <DataArray type="Float64" Name="pe" NumberOfTuples="80" format="appended" RangeMin="12.562471089"         RangeMax="12.644841688"         offset="84"                  />
+    </FieldData>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="Float64" Name="Ca" format="appended" RangeMin="0.0003906208212"      RangeMax="0.0006"               offset="892"                 />
+        <DataArray type="Float64" Name="Cl" format="appended" RangeMin="0.0011983741182"      RangeMax="0.0012"               offset="1344"                />
+        <DataArray type="Float64" Name="H" format="appended" RangeMin="1.0099188654e-07"     RangeMax="1.0117767637e-07"     offset="1760"                />
+        <DataArray type="Float64" Name="K" format="appended" RangeMin="1.0000003266e-12"     RangeMax="0.0003897870855"      offset="2168"                />
+        <DataArray type="Float64" Name="N(5)" format="appended" RangeMin="9.9999999997e-13"     RangeMax="1.6258251135e-06"     offset="2664"                />
+        <DataArray type="Float64" Name="Na" format="appended" RangeMin="1.0000000008e-12"     RangeMax="2.8971274103e-05"     offset="3160"                />
+        <DataArray type="Float64" Name="darcy_velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="3656"                />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="99987.170136"         RangeMax="100000"               offset="3948"                />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="4384"                />
+        <DataArray type="Float64" Name="velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="4444"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="4732"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="5232"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="5396"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="5560"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9A1sTDQTTdPNUw2NUwzMzYHADNCBPA=AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAPAIAAAAAAAA=eF4V0GtIk1EcBvBIp+TykharJW6NLE+mabZmanLmlrXUZGEXzMbKYhKaRfZhySpFM83mSCwXGRTJSoeamF3QNHK5INBVWLhSyM3MYl5Lhhd6zqcfz/Oey/89erc50rqV0Mh7f8QLAkKVknGbW0JoxS5Rv2U7oQ2K+1HHdxO6WFInM8UROtWoHjdQQv2vTxodiYSO5MqaWpMIbRJ/La6SYn1BRESsnNByt8lyWkZol9hScGwPoblHY/o60Ad7jdxKSMb5uqGUN+gdfHO3C5nHn35dCbc8LTt/eS+hnJuq8G1QlDQxNAOPrFW6WuHw85bZjH2ERjdL4gVQSjWf6+FfbZmvFsoTNDW/oNM72qsHcnP0d0UKQjM051Z4wA9Xc9oPwpLpifWx0M05k10In7wPnDoFgziTHXVQ3/etswI+tu989ArW80XmBva9bY3zI3x5O15ngdY0B/kJ/13oddthajz3+xzkhR1qc8GUlelXOPsJ1TiHg5aQ39p6zwYgHxgYK/WBY0tFPutgsVR7aTUUjs7pBFDAVQUHw4UTASEboTV1oH0DzCzyW74ZqneM1obCF18e5IdBj2b3IOtnf3wqZjnNv7qFWd+jymX9tRuZcpbj0v1S2T5fhWfnJmjMvhgrgsm2d1Hsnt/imqwQuOpZYZcQ9o+Gd7P5TIN5BjaPWjmXz+asbDQ+5LH/Utee9IOhMTPegbDZaS/1hhKT+Q4XzldlCZfBEEX5YU9mn65kHu+QmLdkWIT/AXci+AI=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAMQEAAAAAAAA=eF7jMAo10FrpbJ+679uy9pmW9hcDFcF8F4YV+iB60/uleiD6wZ/5OiA6InGnBlj8WY4iiJ65VEwERC97rf9HE0g35PNeAtE7/l7qANEOQSrcINrpxBVTDSD9xiz3vBqQljDdWaUCpBdP612oAKS/PC0IlgLSYW+KBAWAtLH/45r/K5zt70Xy33wCpMvO7XLfBaTdVsnOLgPSFfV1LdxAmsOg/ETScmf7DS5GXyKXOdtPC/U8uXsJUP1vFh2zRc72v2482Oc1z9lebodlkvkMZ/tVW27FXupztjdoMjgX3OBsP6/Rbj1TlrN9j98K5u3uzvZqWWld0sLO9gLJjJKTDznZLznx5vEedyd79ntP6kP7HO27uLScvOod7AMqv8z+o2ZnnzJvKtt8EWt7ANuFihk=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAFgEAAAAAAAA=eF77oRVqoLUy2P4CQ1RmxKJge/fPQWD+j/W+YJpvqhOYDmjTAdOnWpjA9Bv5rfogOs9dC0yz6ZvpgWjTrRN1QLSX9klNEC0QpqIKok0yGmRB9JISCUEQ3ato/FcTSGs8WHkXRC9wtFkDopPmTY8H0dUBv15qAOm7Tj7mIDrFpVZLHUi/83k7UxVIZ62Y6KUMpP16y1/KA+mi7/rM0kDaqnvWD2Eg/fdEzgdOIP1F5eiR3yuC7W89F3jzHEhHaG6MOQOkW+7YZi8D0s1f25YUAelZBhXLdYB0I6ch1/nlwfYbNJXTgoF0DtPeM5uWBds/5na3fr402L7WjyX5/ZJg+zPVyw9cXxxsv4Ynfc55YLgBAJeviQQ=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAEQEAAAAAAAA=eF6bapTa+ECj2u4S+2w+GZlqO+ZrU8F8K5kMMN3l6Q2mj3lKgOn9C9Y3gGiTiyfBtPyDILD4krxrYH7exw/1IPpf/fw6EH352+caEH3t16wqEK3yprEURKdEaeSC6Ll3/eJBtJnNPS8Qva/vvSZYfsbmL/eBdHgG8yIQPcP5uwKIDsq4GHkPSK++f1H1LpDeJdjlfxtIXyiLbLoBpPWOmTlfAdJS6RuvnwXS2xPeah4B0t9Udy3YCqRDm5VuzwfSMmsc/etB/l4avN8PSGucXBHOBqQnCWW1zVIHut97tvEvtWq7+yJLqnmA9Kepb63Xq1TbuUme3r9HqdouwMbqkbZCtZ3X3ubGXNlqOwA2S5WWAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/nQ2kOGZl3E9p9TbyIqLOT96bStl4iN8PRvtVZYeaYo9sqcHhAVXmj0xkc77Ll2qPW2xj5bG6rk9mh2l79fjyD2i0xpzSF3XPbxp8OYrfOU9/iTLlrxm8z1JBygPET4BPhQXTtwpOA4+8+6o8gUnGj4xvvL0emEmPkLA67Xm9TI+hZW+CivVPz4wfq0NwX9KPu63prb141U+Q8BMtc/0YT67BHDee0RtPvhmLCabtXc+VAhMLsoZgz7w1EwxbJ6OPgBTY334bJg+8kxriYRmoz7d/2bd2LGuPnzrw49DMbg+VwuyzbwAwz4atJ8IlMLNPg6p8SSePNc+6tGLxI8Y4j5OMFDIGB7sPukLMg8NzPU+MNL7IjDdAD8msORA8QsKP9S3hrrNERQ/PvNlEiPHHj92Nkwh6EUnPw46+p+W9jA/W71coMa/Nj950ZsuAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/s26K4GZl3E9feD1nd1Guz486QsUWShyPRw/vmAieXM9+OS4Pdc/dj0MwPJLfsN7PTy3hl5UIYM9hmKI7WfNjD1WQtxdeBWXPXhZO2DLKqM9bNGBrmMnsD31NcmfUkm7PXXXp9FY7sY99iIAgxAb0z0oxxpI4ILfPYyGopELtOk9itSbhEK89D3up8yOy4sAPuNuBWNPIAo+nxK62e1pFD69qIxUuZQfPorqd8TEMCg+rDR9KpZaMj5o4Eqk8pc7PjCmpzcgjkQ+9QvjwZtaTj52tFks1TdWPraenTP3H2A+jmqUWbw1Zz4ykHPSAJFwPkehsh/KdHc+T6BAPYl4gD74ktP0UfGGPiqLGGfgsI8+h6njezeylT68rZi4ZmydPtIaI/LyuqM+wv3ivW8Xqj6mM8jAD+uwPlBlz8eOSbU+FrPm4xVquT4rsJ5wAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/u7Fa4GZl3E9lvX4y+pg/j56xGkWcQN0PXClUCniqXk9bEfjmBDSgj1/P31/b8mOPaqE/0LN3Zo9lyS2Y+YYqD2SeeKfB7i1PZukMXpTc8M9ViHafI850T3AMRgvax7ePbeMmwYo/uk9/LyBRjIn9j22/opVq6gCPuuXeOzhFQ8+Hh+XpJmgGT5C9yk8hOwkPiixyUrl7TA+fHuLAn0pOz4Y89fJtZ1FPsZi1pu6ElE+WKUx6nLHWj6nAOAEgdxkPquK6sY4JnA+iik/sAfbeD4K9XmpjwWDPlD6izB/9ow+cMDQrjnxlT7yVVUmJYygPn4KGRun2ag+GBt0pgyVsj4qZfBPQK27PgxK+vvxhcQ+kyp/yUpKzj5aeHHajTfWPha116Y4KOA+CuXF6z8x5z7lzMSjukjwPsAga7qOAfY+aNASYEG6+z4G8p80AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAugAAAAAAAAA=eF4TDdWZpr9kkp1wjByYdvqiBaaX7dIF0x+MjMC0ngCEfvldH0z/y4Dwm1eYQfTNMAfTDKctwHRKtSWYDomD0Nt6IPSeVdZgui3bCkyvkYGoz86Eym+G0GbREHMffDEB012iEHUnTkD076qG2LcuxRBM26TpgekF7hDaoAri/m2MEP8ZXYK4+44NhE6ygMiHnIDQscsgdPRqiP9dfdXAdBqUXrcLQh+KUQHT804qgellcgpgGgDp9YuIAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAJAEAAAAAAAA=eF5jYACCrB8OrNsf7DDK/OGQeNHi3i8gfXhR4Z6vQFq5ZNWsD0C62fVxxWsg/VhMJvwZkHZ5EWL6EEgv3tkrfAdIM3Uf+3gNSCfE/D9/EUjv07VYdwZIy/wv6DkOpKsurMw6BKRvLHzksRdImxZLq+8A0pNdQlg3A+mPor2P1wJp/+dHD64A0mt3/Ju/GEjzdJnXzQPS2dEFMTOB9EmdlVZTgLTav4cS/UC67bzU904g/WxB8NUWIO1e1LO5Hkgvdz46sQpIc4r+yy8F0lnPzPwKgPSp7fk62UBau3MFVxqQ7ol6+CIBSL/VljoeDaT9/gYtDQPSG851NwcCaeEFRxJ9gHRp4V97d5D7nczknIC0jUj+HxsgveDp8lvmQBoAcQ2tMQ==AQAAAAAAAAAAgAAAAAAAAKAAAAAAAAAADAAAAAAAAAA=eF5jYBjcAAAAoAABAQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAtwAAAAAAAAA=eF4TDdWZpr9kkl3+ZE0wffuhIZhmbjEC0xFzIPyqYD0w3VxnBqEvQWjGOAsw/c/cEkyf2Aeh7/2FiPdLW4Hpl5a2YFp5rjmYNpkDkV9kBpG/bgDRN7XVGExrp0DMf34MIr7iFET/GR2IeyZA3fVMXxvizngIn2MqxB/JeRD3MplBxLc764JpXl0IPTsZIi//EuL/6bkaYPqTjgqYvg3lN7KrgmmWTcpg2vKhIpgWipED0wBu+IY5AQAAAAAAAAAAgAAAAAAAANgDAAAAAAAAVgEAAAAAAAA=eF5jYMAOqkXWuT+s2mKPLq7X8fGSb1IChrhWG0i8AEO81TtrT8nkCgzx8k6Q+gYMcda3a4H2tmCIV3qCzOnAFF+uZsixpgdDfHkjyPwJGOKGG/x7p+dNwvTXfZC9UzDEtfl4mLTbp2GI99iA3DMDQ3xB1omyffNnYYgfnQly5xwM8TknWl4H7piHIb6sGuT+BRjiohMVNhRlLMQQ37ES5K9FGOJ5h+uADl2MIX78Jsi/SzDEzX/cUWFrXIoZjyKgcFiGId5iYH3fv3c5hvg2f1D4rMAQf5k3EwhWYprTDQq3VRjiZsu/Bz9euhpDPPMYKDzXYIjPehrKp7tpLYb4P4ZWYDivwxBPlNkCtHg9hvipuaDw34AhnqLb6ZjwFFP87x5QvGzEEJ/us0Ou9TWm+IyboPjahCG+PvX575UfMcU9voLicTOG+MMWsZvnvmOKAwD0N598AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAWAAAAAAAAAA=eF5dxUkCgQAAAMAsoYUUZSkRLf7/QoecZi4TBLMlr3jNIW94yzuOOOaEU97zgTM+cs4Fn/jMJVd84Svf+M41N/zglp/84o7f/OGeBx554i8v/v8AsmsGaA==AQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAWgAAAAAAAAA=eF4txdEGwgAAAMAkk8lMkklmZmaSZGaSyczMJNn/f00P3b3cevW3ceCtQ+8cOfbeBx+d+OSzU2fOXbh05Yuvvvnu2o1bP/x055d7Dx49efbbH3+9+AfO3gZpAQAAAAAAAAAAgAAAAAAAACgAAAAAAAAADAAAAAAAAAA=eF5jZiYOAAAJxAB5
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_400_t_72000.000000.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_400_t_72000.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..ae7c692dbac5fd7e611224f99887850e77bad076
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchangeOutput_ts_400_t_72000.000000.vtu
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+      <DataArray type="Float64" Name="pe" NumberOfTuples="80" format="appended" RangeMin="12.562469116"         RangeMax="12.644840843"         offset="84"                  />
+    </FieldData>
+    <Piece NumberOfPoints="41"                   NumberOfCells="40"                  >
+      <PointData>
+        <DataArray type="Float64" Name="Ca" format="appended" RangeMin="0.00059949413186"     RangeMax="0.0006"               offset="892"                 />
+        <DataArray type="Float64" Name="Cl" format="appended" RangeMin="0.0011999929058"      RangeMax="0.0012"               offset="1288"                />
+        <DataArray type="Float64" Name="H" format="appended" RangeMin="1.0117718638e-07"     RangeMax="1.0117767633e-07"     offset="1640"                />
+        <DataArray type="Float64" Name="K" format="appended" RangeMin="9.9999999999e-13"     RangeMax="9.3841061037e-07"     offset="2004"                />
+        <DataArray type="Float64" Name="N(5)" format="appended" RangeMin="9.9999999997e-13"     RangeMax="7.0949932647e-09"     offset="2500"                />
+        <DataArray type="Float64" Name="Na" format="appended" RangeMin="1e-12"                RangeMax="7.3327670028e-08"     offset="2996"                />
+        <DataArray type="Float64" Name="darcy_velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="3492"                />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="99987.170136"         RangeMax="100000"               offset="3784"                />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="4220"                />
+        <DataArray type="Float64" Name="velocity" format="appended" RangeMin="2.7777777777e-07"     RangeMax="2.7777777788e-07"     offset="4280"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.08"                 offset="4568"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="5068"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="5232"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="5396"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9A1sTDQTTdPNUw2NUwzMzYHADNCBPA=AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAOwIAAAAAAAA=eF4dz2tI02EYBXC3lUplDilnbWkXXXsxmg1CsalvbmtOsqyMZGW12LAIRxqSdBEtJ7Zmhk5zmF2sbW2kiJlUWikhlfkhCi+kzUbJssSyWeYq7fz79OOcD4fn8XOPsp6vI3Q05673dwSh+42S1plYQp1lJxq7JYR2ynOWZSUQqh78RuzxhG7oq+67RAm9Gpln+JhI6KGXutqWZELn9JP6is2Enuu9sTxOTuijQEGzTkboDumQZq+C0GDf2IMO9H86pEXSLYRmhG9M6ELvCrzTNIF8Vt7uNDHuDtAUKgmNMm0NE8MPMbweL7S+8H/XAnUF8e6MFEK/ckUkAnIvJj2xwpGyZF8BZAk+nxmD164oJ5/CmxeqS1arCE3hrPeyYfvjkIad8Ji6ixMHp0sUqlOwtK18QAstrZG2ehjT6HYY4dt0VcVD+HPkepUTejbde/MaxhanHumGlzs1PA8UCLWfhmCS3v/ZNOyXNtRPwC4X7+j8VELtPdOsWWSzTZ3JRd6nbM5ZANOyX82GwYFBq3YJ9OR1lEbAqTwTXwBtIj0/EvZm57etgrfZhTNCOD7Or4qCi2QctQgauZU9a6HbID3I5KZtcWbGrPg2MdPn3uoPZbLPmDbJ7PlZio4zO9Ghh63M7vZy2f01cO7vcP4KeMBVuHAlDA8y7GHuS1Q7cvnwvZDEhkCD2fJlKUw306n/fwT5ZQVDcbFCxYZ1dY6TAXC49rz9F/5mCedJmP8XR5+u+Q69mZW7fsB/a6T2bw==AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAABwEAAAAAAAA=eF47ahhqoLXS2X67ZdKN10uc7X8IQfitx0LAdJ0ihOa5HQSmmeICwHT/Em8wLffdCUzPMjeGqFcWB9Mb1z/VB9GSS/rB9OzK33ogumD9P10Q3SYXrAOiX3S3aoLoQ79uqoBokSx5WRC9oj5TAESf6Fz4WxNIf+q/fhNEq1kaLAbRkWctvEE0n6zfPg0g7ccz8646kJ6g+bRCDUhvm+zgqgJy96aH5xSA9ORnLinSQPresjffhYD0AZkPX1mBtIq4E/enFc72NUqn718G0lPcJM8uA9K5LqnrU4D0vB4NSTYgfdHCvr15ubP9nmX7b91YBvSfz7nZvED627adt92WOtsDAHElh0A=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAA5QAAAAAAAAA=eF77aRFqoLUy2P7qa89JSkA6TgHCj+KD0GF/QsB08CsI7XseQruthtAdjRBayBxC28wIBtMz6oLA9MYVAWC6ucYXTP8J9oDoO2wPphlWGoLpbR+lwPSJ0m/6IPpQ4x4w/Uw2EUw7qV3SA9EHsr7oguiMb/N0QLS/+lQtsHvn1qqD6Lacq0og+ta1e9IgOulWhyCI3rrRlwlEXzZif6MJpFdEbD0Nok2fm88F0b9O/YsC0bs4D/3VAPGTM2tA9I/57WfUgbS+kMsTNSCdaNy1XxXk3+yEHBUg/fNbkZoykAYAvFWEPg==AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAA7wAAAAAAAAA=eF4zYUtrfKBRbcclfffRdSC9o3t3PYif9flSLYhO7Q2pAdH2x19XgeiLZTMqQfSNS/8rQPSEuptgvo1sDFj+wSuOahDdtvACmJZ6+gMsfs9QE8zPk7YBm8eWfhHMd5u5CUwfm50PFl+nUwemm65ngumOXG8wrZ+WA6anZwiCaVuze2BzE8zMwfpXPWMG8wvYp5aBaJk8vkIQrbBwayaIPmf4PQZE7ziW4wXWF/VFF0QHx+mwgGhxP/0j94F0mfCTVBA9U63t7j2QuwOfS4PoQ3vMhe8C6YUuz3ffBvF9RFRvgfxdcWTGDSANAJWsliM=AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/sBlLYGZl3E9HoRlweJ8rz6Oa4+1tZ5xPSUbWChlr3E9the+BRLTcT30grdCJxtyPVTBvHlDp3I9YhLz5+Wvcz02QJDVoZh1PcraEMvpDnk9GE7iY6E9fz1ouIJzexGFPRacEZAbj449KkktLSh2lz0pUY9vKbiiPfvi/Eyyjq49kMFc6D86uT1WHo3bUerEPQjfRga/VtE9eo4W0yCt3D0bVaImwp/nPcAf5YhhYPM9JfDA40mj/z1zC61vRrUJPsL9l48CyxQ+2B4ky7u9ID7qWTuyRNYqPkyfk0YhajU+YkKO9JkCQT5KsokHFeVKPu3uQFMGJ1U+Qq5yCkCKYD4/XVv7Z7FpPu81xDRgy3M+sIR2LTUyfj7QlAeNzb2GPpimUs7g2JA+mpyVI+5qmD6ukqQOjimhPhFZbbLoEKc+eEzJv5rXrD4veZamAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/u6xK4GZl3E9yf5dEAd5Pj41t668y5dxPTSq4ylBmHE9FEL7wDqZcT2CZX8rL5txPfgVTqnynnE9Y1uNKfqlcT28IOPDyrJxPZSFys6wyXE9lrCSWOnxcT04VIKHgDdyPUZn0FpLrnI9CmyQOJl2cz2YnQvnicR0PT57ELRq63Y9jgCFBSpvej1Gm2wKdQ+APf4IG6qQnoQ9Lshgldndiz1y2D6zTaWTPR5os5EhlJw9sEK0pCE4pT2fQZ+Qs+OvPU5+dY4VGbg95nAeEro7wj1bRCBiC4zLPQ49TKeFvNQ9pNSaQnsS3z0ivPIXvybnPY63bXbyJPE9j9HGqVg5+T342MT/UGwCPgoqZQhjswo+8vZl03EtEz4cgiu960IbPlguTWh4ICM+hiZ8YoViKj6wwQrIIMQxPuSRFZC8EDc+I0ZMyGcoPD7shpOSAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAUwEAAAAAAAA=eF4BSAG3/mbqLYGZl3E9xDJqAwqvcz5O366NmphxPWDcbKf0mnE9GAN0X/ifcT2SyhqTEKpxPVgh+dOYvXE9Q+j1j0ricT2K5kcvpCVyPUBCZdzqnnI9AKyoIMl1cz3p7i8AJu10PTj+BSXAdHc9GNI1U43Fez23SufGjoeBPQa0K+zZn4c9dPfXMoLdkD2W6sBl6SyZPfxhavBZXqM9+L7rbqFarj3hIkkiowO4PSXIFOi/EcM9W9I5uiVJzj0G0f2v9f7XPT6yo77W8eI9nDnzS7fI7T2KmY0NDk33PUZ+dTQiIwI+IETTvWkWDD5iuP3zKKEVPoimWv53jyA+9mh/TI4zKT6qkwPBPwszPl0ZpjdYjjw+uRkRGug1RT6Mrp3JICZPPtZY6qysiVY+oNaJusD6Xz4tIPBiMBJmPoIH6d3+O20+L30eLnQYcj7e/5p3AQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAugAAAAAAAAA=eF4TDdWZpr9kkp1wjByYdvqiBaaX7dIF0x+MjMC0ngCEfvldH0z/y4Dwm1eYQfTNMAfTDKctwHRKtSWYDomD0Nt6IPSeVdZgui3bCkyvkYGoz86Eym+G0GbREHMffDEB012iEHUnTkD076qG2LcuxRBM26TpgekF7hDaoAri/m2MEP8ZXYK4+44NhE6ygMiHnIDQscsgdPRqiP9dfdXAdBqUXrcLQh+KUQHT804qgellcgpgGgDp9YuIAQAAAAAAAAAAgAAAAAAAAEgBAAAAAAAAJAEAAAAAAAA=eF5jYACCrB8OrNsf7DDK/OGQeNHi3i8gfXhR4Z6vQFq5ZNWsD0C62fVxxWsg/VhMJvwZkHZ5EWL6EEgv3tkrfAdIM3Uf+3gNSCfE/D9/EUjv07VYdwZIy/wv6DkOpKsurMw6BKRvLHzksRdImxZLq+8A0pNdQlg3A+mPor2P1wJp/+dHD64A0mt3/Ju/GEjzdJnXzQPS2dEFMTOB9EmdlVZTgLTav4cS/UC67bzU904g/WxB8NUWIO1e1LO5Hkgvdz46sQpIc4r+yy8F0lnPzPwKgPSp7fk62UBau3MFVxqQ7ol6+CIBSL/VljoeDaT9/gYtDQPSG851NwcCaeEFRxJ9gHRp4V97d5D7nczknIC0jUj+HxsgveDp8lvmQBoAcQ2tMQ==AQAAAAAAAAAAgAAAAAAAAKAAAAAAAAAADAAAAAAAAAA=eF5jYBjcAAAAoAABAQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAtwAAAAAAAAA=eF4TDdWZpr9kkl3+ZE0wffuhIZhmbjEC0xFzIPyqYD0w3VxnBqEvQWjGOAsw/c/cEkyf2Aeh7/2FiPdLW4Hpl5a2YFp5rjmYNpkDkV9kBpG/bgDRN7XVGExrp0DMf34MIr7iFET/GR2IeyZA3fVMXxvizngIn2MqxB/JeRD3MplBxLc764JpXl0IPTsZIi//EuL/6bkaYPqTjgqYvg3lN7KrgmmWTcpg2vKhIpgWipED0wBu+IY5AQAAAAAAAAAAgAAAAAAAANgDAAAAAAAAVgEAAAAAAAA=eF5jYMAOqkXWuT+s2mKPLq7X8fGSb1IChrhWG0i8AEO81TtrT8nkCgzx8k6Q+gYMcda3a4H2tmCIV3qCzOnAFF+uZsixpgdDfHkjyPwJGOKGG/x7p+dNwvTXfZC9UzDEtfl4mLTbp2GI99iA3DMDQ3xB1omyffNnYYgfnQly5xwM8TknWl4H7piHIb6sGuT+BRjiohMVNhRlLMQQ37ES5K9FGOJ5h+uADl2MIX78Jsi/SzDEzX/cUWFrXIoZjyKgcFiGId5iYH3fv3c5hvg2f1D4rMAQf5k3EwhWYprTDQq3VRjiZsu/Bz9euhpDPPMYKDzXYIjPehrKp7tpLYb4P4ZWYDivwxBPlNkCtHg9hvipuaDw34AhnqLb6ZjwFFP87x5QvGzEEJ/us0Ou9TWm+IyboPjahCG+PvX575UfMcU9voLicTOG+MMWsZvnvmOKAwD0N598AQAAAAAAAAAAgAAAAAAAAIACAAAAAAAAWAAAAAAAAAA=eF5dxUkCgQAAAMAsoYUUZSkRLf7/QoecZi4TBLMlr3jNIW94yzuOOOaEU97zgTM+cs4Fn/jMJVd84Svf+M41N/zglp/84o7f/OGeBx554i8v/v8AsmsGaA==AQAAAAAAAAAAgAAAAAAAAEABAAAAAAAAWgAAAAAAAAA=eF4txdEGwgAAAMAkk8lMkklmZmaSZGaSyczMJNn/f00P3b3cevW3ceCtQ+8cOfbeBx+d+OSzU2fOXbh05Yuvvvnu2o1bP/x055d7Dx49efbbH3+9+AfO3gZpAQAAAAAAAAAAgAAAAAAAACgAAAAAAAAADAAAAAAAAAA=eF5jZiYOAAAJxAB5
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_downstream.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_downstream.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..0f44b032f585f3223cd10bd6dc0e2b814f492dcc
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_downstream.vtu
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="1"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="1"                    RangeMax="1"                    offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="1"                    RangeMax="1"                    offset="24"                  />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.08"                 RangeMax="0.08"                 offset="48"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="92"                  />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="116"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="140"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _CAAAAAAAAAABAAAAAAAAAA==CAAAAAAAAAABAAAAAAAAAA==GAAAAAAAAAB7FK5H4Xq0PwAAAAAAAAAAAAAAAAAAAAA=CAAAAAAAAAAAAAAAAAAAAA==CAAAAAAAAAABAAAAAAAAAA==AQAAAAAAAAAB
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_upstream.vtu b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_upstream.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..03fc506c11d9906bedf6ab450c01822cffca3623
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange_upstream.vtu
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="1"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0"                    RangeMax="0"                    offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="0"                    RangeMax="0"                    offset="24"                  />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0"                    offset="48"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="92"                  />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="116"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="140"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _CAAAAAAAAAAAAAAAAAAAAA==CAAAAAAAAAAAAAAAAAAAAA==GAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA=CAAAAAAAAAAAAAAAAAAAAA==CAAAAAAAAAABAAAAAAAAAA==AQAAAAAAAAAB
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/phreeqc.dat b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/phreeqc.dat
new file mode 100644
index 0000000000000000000000000000000000000000..0873e31a8ef7efd38d7028143691c2d85648c8ed
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CationExchange/phreeqc.dat
@@ -0,0 +1,1837 @@
+# 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  7.81  2.96  -0.46 # supcrt
+# 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  7.81  2.96  -0.46 # supcrt
+H2Sg = HSg- + H+
+	-log_k	-6.994
+	-delta_h 5.30	kcal
+	-analytical  11.17  -0.02386  -3279.0
+	-dw	 1.73e-9
+	-Vm  5.0119  4.9799  3.4765  -2.9849  1.4410 # supcrt
+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
+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	-171.9065	-0.077993	2839.319	71.595
+	-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
+	-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
+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
+
+# Gases from LLNL...
+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.9759
+	-analytic -97.354 -3.1576e-2 1.8285e3 37.44 28.56
+	-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-
+	-analytic -97.354 -3.1576e-2 1.8285e3 37.44 28.56
+	-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
+
+# 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/web/content/docs/benchmarks/reactive-transport/exchange/exchange.md b/web/content/docs/benchmarks/reactive-transport/exchange/exchange.md
new file mode 100644
index 0000000000000000000000000000000000000000..bc53d0b20c67014a7dc9ccf5d6b197d259da4fca
--- /dev/null
+++ b/web/content/docs/benchmarks/reactive-transport/exchange/exchange.md
@@ -0,0 +1,35 @@
++++
+author = "Jaime Garibay-Rodriguez, Renchao Lu, Vanessa Montoya"
+weight = 142
+project = "Parabolic/ComponentTransport/ReactiveTransport/CationExchange/exchange.prj"
+date = "2019-07-08T14:41:09+01:00"
+title = "Transport and Cation Exchange"
+
+[menu]
+  [menu.benchmarks]
+    parent = "Reactive Transport"
+
++++
+
+{{< data-link >}}
+
+## Overview
+
+This benchmark simulates the chemical composition of the effluent from a column containing a cation exchanger (Example 11 in the PHREEQC 3 documentation).
+The following setup is used for the simulation:
+
+{{< img src="../fig1.png" title="Model setup for the simulation of the column with a cation exchanger in OGS6.">}}
+
+Full details of the model setup and parameters are given in the PHREEQC3 example (consulted MAY-2021):
+
+https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqc3-html/phreeqc3-73.htm#50593807_46434
+
+The benchmark uses the `ComponentTransport` process in OGS-6 coupled with the IPhreeqc software (Parkhurst and Appelo,2013). The results show good agreement between codes. More details about the implementation of the `ComponentTransport` process in OGS-6 can be found in  [HC-Process.pdf](/docs/benchmarks/hydro-component/HC-Process.pdf).
+
+{{< img src="../fig2.png" title="Comparison between PHREEQC and OGS6 of simulated concentrations of solutes at time = 18,000 s reacting with an exchanger.">}}
+
+{{< data-link >}}
+
+## References
+
+Parkhurst, D.L., Appelo, C.A.J., 2013. Description of Input and Examples for PHREEQC Version 3 - a Computer Program for Speciation, Batch-reaction, One-dimensional Transport, and Inverse Geochemical Calculations.
\ No newline at end of file
diff --git a/web/content/docs/benchmarks/reactive-transport/exchange/fig1.png b/web/content/docs/benchmarks/reactive-transport/exchange/fig1.png
new file mode 100644
index 0000000000000000000000000000000000000000..84b3c2553ac4b6539b0b30008fd19531e9262221
Binary files /dev/null and b/web/content/docs/benchmarks/reactive-transport/exchange/fig1.png differ
diff --git a/web/content/docs/benchmarks/reactive-transport/exchange/fig2.png b/web/content/docs/benchmarks/reactive-transport/exchange/fig2.png
new file mode 100644
index 0000000000000000000000000000000000000000..12cca9022ce5b6f95107a6173121f69435ab99a6
Binary files /dev/null and b/web/content/docs/benchmarks/reactive-transport/exchange/fig2.png differ