diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index 998581e63582affc2b027484444d60fae45a9c17..83f3afaa1b499129a0e4e94c6865286c3ef820d3 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -104,9 +104,7 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
         config.getConfigSubtreeOptional("rates"));
 
     // surface
-    auto surface = PhreeqcIOData::createSurface(
-        //! \ogs_file_param{prj__chemical_system__surface}
-        config.getConfigSubtreeOptional("surface"));
+    auto const& surface = chemical_system->surface;
 
     // exchange
     auto const& exchangers = chemical_system->exchangers;
@@ -147,8 +145,8 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
     return std::make_unique<PhreeqcIOData::PhreeqcIO>(
         *linear_solver, std::move(project_file_name),
         std::move(path_to_database), std::move(chemical_system),
-        std::move(reaction_rates), std::move(surface), std::move(user_punch),
-        std::move(output), std::move(dump), std::move(knobs));
+        std::move(reaction_rates), std::move(user_punch), std::move(output),
+        std::move(dump), std::move(knobs));
 }
 
 template <>
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index c2e1e1e138714dd5cd7dde6bfed080203a3dccb9..f754dd71f42c2261a2e04e07b55757f4cea9646d 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -42,6 +42,14 @@ namespace PhreeqcIOData
 {
 namespace
 {
+template <class... Ts>
+struct overloaded : Ts...
+{
+    using Ts::operator()...;
+};
+template <class... Ts>
+overloaded(Ts...)->overloaded<Ts...>;
+
 template <typename DataBlock>
 std::ostream& operator<<(std::ostream& os,
                          std::vector<DataBlock> const& data_blocks)
@@ -164,20 +172,20 @@ void setReactantMolality(Reactant& reactant,
         (*reactant.molality)[chemical_system_id];
 }
 
-template <typename Exchanger>
-void initializeExchangerMolality(Exchanger& exchanger,
-                                 GlobalIndexType const& chemical_system_id,
-                                 MaterialPropertyLib::Phase const& solid_phase,
-                                 ParameterLib::SpatialPosition const& pos,
-                                 double const t)
+template <typename Site>
+void initializeSiteMolality(Site& site,
+                            GlobalIndexType const& chemical_system_id,
+                            MaterialPropertyLib::Phase const& solid_phase,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t)
 {
-    auto const& solid_constituent = solid_phase.component(exchanger.name);
+    auto const& solid_constituent = solid_phase.component(site.name);
 
     auto const molality =
         solid_constituent[MaterialPropertyLib::PropertyType::molality]
             .template initialValue<double>(pos, t);
 
-    (*exchanger.molality)[chemical_system_id] = molality;
+    (*site.molality)[chemical_system_id] = molality;
 }
 
 template <typename Reactant>
@@ -253,7 +261,6 @@ PhreeqcIO::PhreeqcIO(GlobalLinearSolver& linear_solver,
                      std::string&& database,
                      std::unique_ptr<ChemicalSystem>&& chemical_system,
                      std::vector<ReactionRate>&& reaction_rates,
-                     std::vector<SurfaceSite>&& surface,
                      std::unique_ptr<UserPunch>&& user_punch,
                      std::unique_ptr<Output>&& output,
                      std::unique_ptr<Dump>&& dump,
@@ -263,7 +270,6 @@ PhreeqcIO::PhreeqcIO(GlobalLinearSolver& linear_solver,
       _database(std::move(database)),
       _chemical_system(std::move(chemical_system)),
       _reaction_rates(std::move(reaction_rates)),
-      _surface(std::move(surface)),
       _user_punch(std::move(user_punch)),
       _output(std::move(output)),
       _dump(std::move(dump)),
@@ -341,8 +347,18 @@ void PhreeqcIO::initializeChemicalSystemConcrete(
 
     for (auto& exchanger : _chemical_system->exchangers)
     {
-        initializeExchangerMolality(exchanger, chemical_system_id, solid_phase,
-                                    pos, t);
+        initializeSiteMolality(exchanger, chemical_system_id, solid_phase, pos,
+                               t);
+    }
+
+    for (auto& surface_site : _chemical_system->surface)
+    {
+        if (auto const surface_site_ptr =
+                std::get_if<MoleBasedSurfaceSite>(&surface_site))
+        {
+            initializeSiteMolality(*surface_site_ptr, chemical_system_id,
+                                   solid_phase, pos, t);
+        }
     }
 }
 
@@ -527,7 +543,7 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
                << "\n";
         }
 
-        auto const& surface = phreeqc_io._surface;
+        auto const& surface = phreeqc_io._chemical_system->surface;
         if (!surface.empty())
         {
             os << "SURFACE " << chemical_system_id + 1 << "\n";
@@ -536,9 +552,44 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
                     ? chemical_system_id + 1
                     : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
             os << "-equilibrate with solution " << aqueous_solution_id << "\n";
-            os << "-sites_units DENSITY"
-               << "\n";
-            os << surface << "\n";
+
+            // print unit
+            if (std::holds_alternative<DensityBasedSurfaceSite>(
+                    surface.front()))
+            {
+                os << "-sites_units density"
+                   << "\n";
+            }
+            else
+            {
+                os << "-sites_units absolute"
+                   << "\n";
+            }
+
+            for (auto const& surface_site : surface)
+            {
+                std::visit(
+                    overloaded{[&os](DensityBasedSurfaceSite const& s) {
+                                   os << s.name << " " << s.site_density << " "
+                                      << s.specific_surface_area << " "
+                                      << s.mass << "\n";
+                               },
+                               [&os, chemical_system_id](
+                                   MoleBasedSurfaceSite const& s) {
+                                   os << s.name << " "
+                                      << (*s.molality)[chemical_system_id]
+                                      << "\n";
+                               }},
+                    surface_site);
+            }
+
+            // overlook the effect of the buildup of charges onto the surface
+            if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
+            {
+                os << "-no_edl"
+                   << "\n";
+            }
+
             os << "SAVE solution " << chemical_system_id + 1 << "\n";
         }
 
@@ -617,7 +668,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
     auto const& output = *phreeqc_io._output;
     auto const& dropped_item_ids = output.dropped_item_ids;
 
-    auto const& surface = phreeqc_io._surface;
+    auto const& surface = phreeqc_io._chemical_system->surface;
     auto const& exchangers = phreeqc_io._chemical_system->exchangers;
 
     if (!surface.empty() && !exchangers.empty())
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index bca62ca832b530ebd32401b5eaf6575b45e66f47..5b48e35d61d3f6391ae74b5fc46aaa3838a83a33 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -27,7 +27,6 @@ namespace PhreeqcIOData
 struct ChemicalSystem;
 struct ReactionRate;
 struct Output;
-struct SurfaceSite;
 struct Dump;
 struct UserPunch;
 
@@ -39,7 +38,6 @@ public:
               std::string&& database,
               std::unique_ptr<ChemicalSystem>&& chemical_system,
               std::vector<ReactionRate>&& reaction_rates,
-              std::vector<SurfaceSite>&& surface,
               std::unique_ptr<UserPunch>&& user_punch,
               std::unique_ptr<Output>&& output,
               std::unique_ptr<Dump>&& dump,
@@ -111,7 +109,6 @@ private:
     std::string const _database;
     std::unique_ptr<ChemicalSystem> _chemical_system;
     std::vector<ReactionRate> const _reaction_rates;
-    std::vector<SurfaceSite> const _surface;
     std::unique_ptr<UserPunch> _user_punch;
     std::unique_ptr<Output> const _output;
     std::unique_ptr<Dump> const _dump;
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
index c3b704f0ee0439aa9f94147fc282fa1c9badec67..b1c4956410fba618679e4531e8ed4ed8e4075e80 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
@@ -56,6 +56,15 @@ void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
     {
         exchanger.molality->resize(num_chemical_systems);
     }
+
+    for (auto& surface_site : surface)
+    {
+        if (auto const surface_site_ptr =
+                std::get_if<MoleBasedSurfaceSite>(&surface_site))
+        {
+            surface_site_ptr->molality->resize(num_chemical_systems);
+        }
+    }
 }
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.h b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
index 9e50e9cb5ff2ba2160f78cb412f2a564f5a1fb57..22b209ca175badf3ddde1d84126dbb055b6eb19c 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
@@ -11,12 +11,14 @@
 #pragma once
 
 #include <memory>
+#include <variant>
 #include <vector>
 
 #include "AqueousSolution.h"
 #include "EquilibriumReactant.h"
 #include "Exchange.h"
 #include "KineticReactant.h"
+#include "Surface.h"
 
 namespace ChemistryLib
 {
@@ -27,11 +29,14 @@ struct ChemicalSystem
     ChemicalSystem(std::unique_ptr<AqueousSolution>&& aqueous_solution_,
                    std::vector<KineticReactant>&& kinetic_reactants_,
                    std::vector<EquilibriumReactant>&& equilibrium_reactants_,
-                   std::vector<ExchangeSite>&& exchangers_)
+                   std::vector<ExchangeSite>&& exchangers_,
+                   std::vector<std::variant<DensityBasedSurfaceSite,
+                                            MoleBasedSurfaceSite>>&& surface_)
         : aqueous_solution(std::move(aqueous_solution_)),
           kinetic_reactants(std::move(kinetic_reactants_)),
           equilibrium_reactants(std::move(equilibrium_reactants_)),
-          exchangers(std::move(exchangers_))
+          exchangers(std::move(exchangers_)),
+          surface(std::move(surface_))
     {
     }
 
@@ -41,6 +46,8 @@ struct ChemicalSystem
     std::vector<KineticReactant> kinetic_reactants;
     std::vector<EquilibriumReactant> equilibrium_reactants;
     std::vector<ExchangeSite> exchangers;
+    std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
+        surface;
 };
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp
index 6319c899a41d0e290c14e2449a75f34e91618911..1aef5665d1ea563ab0d045d99539690318fdc29e 100644
--- a/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp
@@ -16,9 +16,11 @@
 #include "CreateEquilibriumReactants.h"
 #include "CreateExchange.h"
 #include "CreateKineticReactant.h"
+#include "CreateSurface.h"
 #include "EquilibriumReactant.h"
 #include "Exchange.h"
 #include "KineticReactant.h"
+#include "Surface.h"
 
 namespace ChemistryLib
 {
@@ -47,10 +49,16 @@ std::unique_ptr<ChemicalSystem> createChemicalSystem(
         //! \ogs_file_param{prj__chemical_system__exchangers}
         config.getConfigSubtreeOptional("exchangers"), mesh);
 
+    // surface
+    auto surface = createSurface(
+        //! \ogs_file_param{prj__chemical_system__surface}
+        config.getConfigSubtreeOptional("surface"), mesh);
+
     return std::make_unique<ChemicalSystem>(std::move(aqueous_solution),
                                             std::move(kinetic_reactants),
                                             std::move(equilibrium_reactants),
-                                            std::move(exchangers));
+                                            std::move(exchangers),
+                                            std::move(surface));
 }
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateSurface.cpp b/ChemistryLib/PhreeqcIOData/CreateSurface.cpp
index 9c526bac265bb23074c753cebed6d7b59505e4fa..161767afb0b0261333c0fc31cf288d5cc7d74cfe 100644
--- a/ChemistryLib/PhreeqcIOData/CreateSurface.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateSurface.cpp
@@ -8,48 +8,79 @@
  *
  */
 
-#include <optional>
+#include "CreateSurface.h"
 
 #include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
 #include "Surface.h"
 
 namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-std::vector<SurfaceSite> createSurface(
-    std::optional<BaseLib::ConfigTree> const& config)
+std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
+createSurface(std::optional<BaseLib::ConfigTree> const& config,
+              MeshLib::Mesh& mesh)
 {
     if (!config)
     {
         return {};
     }
 
-    std::vector<SurfaceSite> surface;
-    for (auto const& site_config :
-         //! \ogs_file_param{prj__chemical_system__surface__site}
-         config->getConfigSubtreeList("site"))
+    std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
+        surface;
+
+    auto const surface_site_unit =
+        //! \ogs_file_attr{prj__chemical_system__surface__site_unit}
+        config->getConfigAttribute<std::string>("site_unit", "mole");
+
+    if (surface_site_unit == "density")
     {
-        //! \ogs_file_param{prj__chemical_system__surface__site__name}
-        auto name = site_config.getConfigParameter<std::string>("name");
+        for (auto const& site_config :
+             //! \ogs_file_param{prj__chemical_system__surface__site}
+             config->getConfigSubtreeList("site"))
+        {
+            //! \ogs_file_param{prj__chemical_system__surface__site__name}
+            auto name = site_config.getConfigParameter<std::string>("name");
+
+            auto const site_density =
+                //! \ogs_file_param{prj__chemical_system__surface__site__site_density}
+                site_config.getConfigParameter<double>("site_density");
+
+            auto const specific_surface_area =
+                //! \ogs_file_param{prj__chemical_system__surface__site__specific_surface_area}
+                site_config.getConfigParameter<double>("specific_surface_area");
+
+            auto const mass =
+                //! \ogs_file_param{prj__chemical_system__surface__site__mass}
+                site_config.getConfigParameter<double>("mass");
 
-        auto const site_density =
-            //! \ogs_file_param{prj__chemical_system__surface__site__site_density}
-            site_config.getConfigParameter<double>("site_density");
+            surface.push_back(DensityBasedSurfaceSite(
+                std::move(name), site_density, specific_surface_area, mass));
+        }
+
+        return surface;
+    }
+
+    if (surface_site_unit == "mole")
+    {
+        for (auto const& site_config :
+             //! \ogs_file_param{prj__chemical_system__surface__site}
+             config->getConfigSubtreeList("site"))
+        {
+            //! \ogs_file_param{prj__chemical_system__surface__site__name}
+            auto name = site_config.getConfigParameter<std::string>("name");
 
-        auto const specific_surface_area =
-            //! \ogs_file_param{prj__chemical_system__surface__site__specific_surface_area}
-            site_config.getConfigParameter<double>("specific_surface_area");
+            auto const molality = MeshLib::getOrCreateMeshProperty<double>(
+                mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
-        auto const mass =
-            //! \ogs_file_param{prj__chemical_system__surface__site__mass}
-            site_config.getConfigParameter<double>("mass");
+            surface.push_back(MoleBasedSurfaceSite(std::move(name), molality));
+        }
 
-        surface.emplace_back(
-            std::move(name), site_density, specific_surface_area, mass);
+        return surface;
     }
 
-    return surface;
+    OGS_FATAL("Surface site unit should be either of 'density' or 'mole'.");
 }
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateSurface.h b/ChemistryLib/PhreeqcIOData/CreateSurface.h
index 849dd6fe72ffa2fcd1a2c1e6c1fed7184220f931..392b435cc5d409f77f45cc49808b2544597e0250 100644
--- a/ChemistryLib/PhreeqcIOData/CreateSurface.h
+++ b/ChemistryLib/PhreeqcIOData/CreateSurface.h
@@ -11,6 +11,7 @@
 #pragma once
 
 #include <optional>
+#include <variant>
 #include <vector>
 
 namespace BaseLib
@@ -18,13 +19,20 @@ namespace BaseLib
 class ConfigTree;
 }
 
+namespace MeshLib
+{
+class Mesh;
+}
+
 namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-struct SurfaceSite;
+struct DensityBasedSurfaceSite;
+struct MoleBasedSurfaceSite;
 
-std::vector<SurfaceSite> createSurface(
-    std::optional<BaseLib::ConfigTree> const& config);
+std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
+createSurface(std::optional<BaseLib::ConfigTree> const& config,
+              MeshLib::Mesh& mesh);
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/Surface.cpp b/ChemistryLib/PhreeqcIOData/Surface.cpp
deleted file mode 100644
index 53fbafc54beebdd4be65db7bfc49b74eecd8da74..0000000000000000000000000000000000000000
--- a/ChemistryLib/PhreeqcIOData/Surface.cpp
+++ /dev/null
@@ -1,28 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "Surface.h"
-
-#include <ostream>
-
-namespace ChemistryLib
-{
-namespace PhreeqcIOData
-{
-std::ostream& operator<<(std::ostream& os, SurfaceSite const& surface_site)
-{
-    os << surface_site.name << " " << surface_site.site_density << " "
-       << surface_site.specific_surface_area << " " << surface_site.mass
-       << "\n";
-
-    return os;
-}
-}  // namespace PhreeqcIOData
-}  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/Surface.h b/ChemistryLib/PhreeqcIOData/Surface.h
index 4edb777a0c562f8ee0a7f630b3d67dc4977b5e3b..aae9ae54ed800bfc2aff569c3b7e8ebb4fd9d079 100644
--- a/ChemistryLib/PhreeqcIOData/Surface.h
+++ b/ChemistryLib/PhreeqcIOData/Surface.h
@@ -10,24 +10,20 @@
 
 #pragma once
 
-#include <iosfwd>
 #include <string>
 
-namespace BaseLib
-{
-class ConfigTree;
-}
+#include "MeshLib/PropertyVector.h"
 
 namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-struct SurfaceSite
+struct DensityBasedSurfaceSite
 {
-    SurfaceSite(std::string name_,
-                double const site_density_,
-                double const specific_surface_area_,
-                double const mass_)
+    DensityBasedSurfaceSite(std::string name_,
+                            double const site_density_,
+                            double const specific_surface_area_,
+                            double const mass_)
         : name(std::move(name_)),
           site_density(site_density_),
           specific_surface_area(specific_surface_area_),
@@ -35,13 +31,22 @@ struct SurfaceSite
     {
     }
 
-    friend std::ostream& operator<<(std::ostream& os,
-                                    SurfaceSite const& surface_site);
-
     std::string const name;
     double const site_density;
     double const specific_surface_area;
     double const mass;
 };
+
+struct MoleBasedSurfaceSite
+{
+    MoleBasedSurfaceSite(std::string name_,
+                         MeshLib::PropertyVector<double>* const molality_)
+        : name(std::move(name_)), molality(molality_)
+    {
+    }
+
+    std::string const name;
+    MeshLib::PropertyVector<double>* molality;
+};
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/Documentation/ProjectFile/prj/chemical_system/surface/a_site_unit.md b/Documentation/ProjectFile/prj/chemical_system/surface/a_site_unit.md
new file mode 100644
index 0000000000000000000000000000000000000000..5bcba7ad2acff716d2a3c72cec56c589e53eb3c4
--- /dev/null
+++ b/Documentation/ProjectFile/prj/chemical_system/surface/a_site_unit.md
@@ -0,0 +1 @@
+Define site unit for surface complexation model. Either `density` to use sites per square nanometre or `mole` to use moles per kilogram of water are allowed.
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index 9317c623a967b94fc14352ef8517a374feb02310..281e5468f586258926d0b1cc49d785ccf0670067 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -816,6 +816,7 @@ if (NOT OGS_USE_MPI)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/KineticReactant/1d_isofrac_flag_formula.prj RUNTIME 20)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/KineticReactant_AllAsComponents/KineticReactant2.prj RUNTIME 60)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj RUNTIME 60)
+    OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption_SitesMoles.prj RUNTIME 33)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/LookupTable/RadionuclideSorption.prj RUNTIME 10)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/EquilibriumPhase/calciteDissolvePrecipitateOnly.prj RUNTIME 25)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/EquilibriumPhase/calcitePorosityChange.prj RUNTIME 25)
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj
index 82988e7ddba4a5e40c385663ff8859ab74d8c516..b7a295da4cd7ae3f9c39f4a55c5a2b9f0b854bbf 100644
--- a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption.prj
@@ -464,7 +464,7 @@
                 <saturation_index>0.0</saturation_index>
             </phase_component>
         </equilibrium_reactants>
-        <surface>
+        <surface site_unit="density">
             <site>
                 <!--Qz  = Mineral phase: Quartz, amorphous Silica-->
 		        <name>QzOH</name>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption_SitesMoles.prj b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption_SitesMoles.prj
new file mode 100644
index 0000000000000000000000000000000000000000..349c89c911634c142aa1bea93fa368b44fdf58c2
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/RadionuclideSorption_SitesMoles.prj
@@ -0,0 +1,859 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>RadionuclideSorption.vtu</mesh>
+        <mesh>RadionuclideSorption_ReactiveDomain.vtu</mesh>
+        <mesh>RadionuclideSorption_left.vtu</mesh>
+        <mesh>RadionuclideSorption_right.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>hc</name>
+            <type>ComponentTransport</type>
+            <integration_order>2</integration_order>
+            <coupling_scheme>staggered</coupling_scheme>
+            <process_variables>
+                <concentration>Na</concentration>
+                <concentration>Cl</concentration>
+                <concentration>Ca</concentration>
+                <concentration>C(4)</concentration>
+                <concentration>H</concentration>
+                <!--Radionuclide-->
+                <concentration>U(6)</concentration>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <specific_body_force>0 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>Na</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>Cl</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>Ca</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>C(4)</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>H</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>U(6)</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                    <properties>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value>1e3</value>
+                        </property>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value>1e-3</value>
+                        </property>
+                    </properties>
+                </phase>
+                <phase>
+                    <type>Solid</type>
+                    <components>
+                        <component>
+                            <name>Calcite</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>10</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>AlbOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>8.538e-04</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>GibbOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>2.236e-05</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>HemOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>5.285e-05</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>KaoOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>4.269e-05</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>MusOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>3.497e-04</value>
+                                </property>
+                            </properties>
+                        </component>
+                        <component>
+                            <name>QzOH</name>
+                            <properties>
+                                <property>
+                                    <name>molality</name>
+                                    <type>Constant</type>
+                                    <value>2.419e-04</value>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>permeability</name>
+                    <type>Parameter</type>
+                    <parameter_name>kappa</parameter_name>
+                </property>
+                <property>
+                    <name>porosity</name>
+                    <type>Parameter</type>
+                    <parameter_name>porosity</parameter_name>
+                </property>
+                <property>
+                    <name>longitudinal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0.2</value>
+                </property>
+                <property>
+                    <name>transversal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0.2</value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <time_loop>
+        <global_process_coupling>
+            <max_iter>6</max_iter>
+            <convergence_criteria>
+                <!-- convergence criterion for the first process (p) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (Na) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (Cl) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (Ca) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (C(4)) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (H) -->
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-14</reltol>
+                </convergence_criterion>
+                <!-- convergence criterion for the second process (U(6)) -->
+                <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 hydraulic equation -->
+             <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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (C(4)) -->
+             <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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+           <!-- convergence criterion for component transport equation (U(6)) -->
+             <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>115000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>115</repeat>
+                            <delta_t>1000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>RadionuclideSorption</prefix>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+            <timesteps>
+                <pair>
+                    <repeat>5</repeat>
+                    <each_steps>23</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>Na</variable>
+                <variable>Cl</variable>
+                <variable>Ca</variable>
+                <variable>C(4)</variable>
+                <variable>H</variable>
+                <variable>U(6)</variable>
+                <variable>pressure</variable>
+                <variable>darcy_velocity</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <chemical_system chemical_solver="Phreeqc">
+        <mesh>RadionuclideSorption_ReactiveDomain</mesh>
+        <linear_solver>general_linear_solver</linear_solver>
+        <database>PSINA_12_07_110615_DAV_s_hzdr.dat</database>
+        <solution>
+            <temperature>14</temperature>
+            <pressure>1</pressure>
+            <pe>4</pe>
+            <components>
+                <component>Na</component>
+                <component>Cl</component>
+                <component>Ca</component>
+                <component>C(4)</component>
+                <!--Radionuclides-->
+                <component>U(6)</component>
+            </components>
+        </solution>
+        <equilibrium_reactants>
+            <phase_component>
+		    <name>Calcite</name>
+                <saturation_index>0.0</saturation_index>
+            </phase_component>
+        </equilibrium_reactants>
+        <surface site_unit="mole">
+            <site>
+                <!--Qz  = Mineral phase: Quartz, amorphous Silica-->
+		        <name>QzOH</name>
+            </site>
+            <site>
+                <!--Alb = Mineral group: Feldspar (albite)-->
+		        <name>AlbOH</name>
+            </site>
+            <site>
+                <!--Mus  = Mineral group: Mica (muscovite)-->
+		        <name>MusOH</name>
+            </site>
+            <site>
+                <!--Hem  = Mineral group: Fe(III)-oxids/-hydroxids (hematite)-->
+		        <name>HemOH</name>
+            </site>
+            <site>
+                <!--Gibb  = Mineral group: Al-hydroxids (gibbsite)-->
+		        <name>GibbOH</name>
+            </site>
+            <site>
+                <!--Kao = Mineral group: 2-layer-clay minerals (kaolinite)-->
+		        <name>KaoOH</name>
+            </site>
+        </surface>
+        <user_punch>
+	        <headline>Kd_U(6)</headline>
+            <statement>k1=(mol("QzOUO2CO3-") + mol("QzOUO2+")+ mol("QzOUO2(OH)"))</statement>
+            <statement>k2=(mol("AlbOUO2+"))</statement>
+            <statement>k3=(mol("MusOUO2+") + mol("MusOUO2CO3-"))</statement>
+            <statement>k4=(mol("HemOUO2CO3-") + mol("HemOUO2+"))</statement>
+            <statement>k5=(mol("(GibbO)2UO2") + mol("GibbOUO2(OH)")+ mol("GibbO(UO2)3(OH)5") + mol("GibbOUO2+"))</statement>
+            <statement>k6=(mol("KaoOUO2+") + mol("KaoOUO2(OH)") + mol("KaoOUO2CO3-"))</statement>
+            <!-- -->
+            <statement>g_Tot=10388   # Kd = (mol_sorb/g_tot)/(mol_aq)  -->  Kd in l/g = m³/kg</statement>
+            <statement>if TOT("U(6)") &gt; 0 then kd_U(6) = ((k1+k2+k3+k4+k5+k6)/TOT("U(6)")/g_tot) else kd_U(6) = 0 </statement>
+            <statement>PUNCH kd_U(6)</statement>
+        </user_punch>
+        <knobs>
+            <max_iter>250</max_iter>
+            <relative_convergence_tolerance>1e-6</relative_convergence_tolerance>
+            <tolerance>1e-20</tolerance>
+            <step_size>5</step_size>
+            <scaling>1</scaling>
+        </knobs>
+    </chemical_system>
+    <parameters>
+        <parameter>
+            <name>kappa</name>
+            <type>Constant</type>
+            <values>1.019368e-11</values>
+        </parameter>
+        <parameter>
+            <name>porosity</name>
+            <type>Constant</type>
+            <value>0.2</value>
+        </parameter>
+        <parameter>
+            <name>retardation</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>decay</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p_Neumann_left</name>
+            <type>Constant</type>
+            <value>3e-2</value>
+        </parameter>
+        <parameter>
+            <name>p_Dirichlet_right</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>c0_Na</name>
+            <type>Constant</type>
+            <value>1e-3</value>
+        </parameter>
+        <parameter>
+            <name>c0_Cl</name>
+            <type>Constant</type>
+            <value>6e-2</value>
+        </parameter>
+        <parameter>
+            <name>c0_Ca</name>
+            <type>Constant</type>
+            <value>2.502e-2</value>
+        </parameter>
+        <parameter>
+            <name>c0_C(4)</name>
+            <type>Constant</type>
+            <value>1e-3</value>
+        </parameter>
+        <parameter>
+            <name>c0_H</name>
+            <type>Constant</type>
+            <!--pH=7-->
+            <value>1e-7</value>
+        </parameter>
+        <parameter>
+            <name>c0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>c_left_Na</name>
+            <type>Constant</type>
+            <value>1e-3</value>
+        </parameter>
+        <parameter>
+            <name>c_left_Cl</name>
+            <type>Constant</type>
+            <value>6e-2</value>
+        </parameter>
+        <parameter>
+            <name>c_left_Ca</name>
+            <type>Constant</type>
+            <value>2.502e-2</value>
+        </parameter>
+        <parameter>
+            <name>c_left_C(4)</name>
+            <type>Constant</type>
+            <value>1e-3</value>
+        </parameter>
+        <parameter>
+            <name>c_left_H</name>
+            <type>Constant</type>
+            <!--pH=7-->
+            <value>1e-7</value>
+        </parameter>
+        <parameter>
+            <name>c_spatial</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>c_left_default</name>
+            <type>CurveScaled</type>
+            <curve>c_temporal_default</curve>
+            <parameter>c_spatial</parameter>
+        </parameter>
+    </parameters>
+    <curves>
+        <curve>
+            <name>c_temporal_default</name>
+            <coords>0 10000 10000.001 115000</coords>
+            <values>1e-8 1e-8 0 0</values>
+        </curve>
+    </curves>
+    <process_variables>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>RadionuclideSorption_left</mesh>
+                    <type>Neumann</type>
+                    <parameter>p_Neumann_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>RadionuclideSorption_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet_right</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>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_Na</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>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_Cl</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>Ca</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_Ca</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_Ca</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>C(4)</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0_C(4)</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_C(4)</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>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_H</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <!--Radionuclides-->
+        <process_variable>
+            <name>U(6)</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>RadionuclideSorption_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_left_default</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>10</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p ilut -tol 1e-8 -maxiter 20000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-12</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>hc</prefix>
+                <parameters>-hc_ksp_type bcgs -hc_pc_type bjacobi -hc_ksp_rtol 1e-8 -hc_ksp_max_it 20000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <regex>RadionuclideSorption_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>RadionuclideSorption_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>darcy_velocity</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-6</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>RadionuclideSorption_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>
+        <vtkdiff>
+            <regex>RadionuclideSorption_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>RadionuclideSorption_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>RadionuclideSorption_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>C(4)</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <regex>RadionuclideSorption_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>RadionuclideSorption_ts_[0-9]*_t_[0-9]*.000000.vtu</regex>
+            <field>U(6)</field>
+            <absolute_tolerance>1e-10</absolute_tolerance>
+            <relative_tolerance>1e-16</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>