diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 60954b084caadd2969f0f314b7e9bba1f74c22de..a038048709e3d97481938d59a05ca7d014e48216 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -34,8 +34,8 @@
 
 #include "ProcessLib/UncoupledProcessesTimeLoop.h"
 
-#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
 #include "ProcessLib/ComponentTransport/CreateComponentTransportProcess.h"
+#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
 #include "ProcessLib/HT/CreateHTProcess.h"
 #include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
 #include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
@@ -43,6 +43,7 @@
 #include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
 #include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
 #include "ProcessLib/PhaseField/CreatePhaseFieldProcess.h"
+#include "ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h"
 #include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
@@ -420,6 +421,14 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                     break;
             }
         }
+        else if (type == "RichardsComponentTransport")
+        {
+            process = ProcessLib::RichardsComponentTransport::
+                createRichardsComponentTransportProcess(
+                    *_mesh_vec[0], std::move(jacobian_assembler),
+                    _process_variables, _parameters, integration_order,
+                    process_config);
+        }
         else if (type == "SMALL_DEFORMATION")
         {
             switch (_mesh_vec[0]->getDimension())
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/c_RichardsComponentTransport.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/c_RichardsComponentTransport.md
new file mode 100644
index 0000000000000000000000000000000000000000..a3a79b32e1a0007f99fb50c1c0487e5660733add
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/c_RichardsComponentTransport.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::RichardsComponentTransport::RichardsComponentTransportProcess
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..d7849410e30e3af70dde2746445407b46f6c9d61
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/i_porous_medium.md
@@ -0,0 +1,3 @@
+The id of the porous medium as an integer between 0 and the number of material
+groups minus one. There must be a corresponding data array 'MaterialIDs' in the
+mesh file.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/a_id.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/a_id.md
new file mode 100644
index 0000000000000000000000000000000000000000..c1dae9f80b40a7c5ba5eb2fc7fae385bf250a298
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/a_id.md
@@ -0,0 +1,2 @@
+It specifies the material ID of a porous medium. The material IDs of elements
+must be given in the mesh's cell data array "MaterialIDs".
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..bb330dc056ceff8400e948dfb52a53fb0f44252a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/i_porous_medium.md
@@ -0,0 +1 @@
+A tag for the properties of a porous medium with material ID.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_capillary_pressure.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_capillary_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..19d854870d0c3537197cc8edb349b03503fcbede
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_capillary_pressure.md
@@ -0,0 +1 @@
+Capillary pressure model.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_permeability.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..075c8a48fde6b9a1d7682d9cc08e6bcaa9feeb93
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_permeability.md
@@ -0,0 +1 @@
+A tag for the relative permeability model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_porosity.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_porosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..d90fcba994c498eb7f7c82b6052cf4619ef7735b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_porosity.md
@@ -0,0 +1 @@
+A tag for the porosity model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_relative_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..fd2a51343fb9e4a7a8ce38bdc88c31f400c1b7e7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_relative_permeability.md
@@ -0,0 +1 @@
+A tag for the relative permeability model.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_storage.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_storage.md
new file mode 100644
index 0000000000000000000000000000000000000000..756a75cbc74b625cca84e3b6b15dcc29b36026b5
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/porous_medium/porous_medium/t_storage.md
@@ -0,0 +1 @@
+A tag for the storage model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/i_process_variables.md
new file mode 100644
index 0000000000000000000000000000000000000000..322cded9564f725b17fb8d83a1991c75376ea8a6
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/i_process_variables.md
@@ -0,0 +1 @@
+The process variable for pressure and concentration.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_concentration.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_concentration.md
new file mode 100644
index 0000000000000000000000000000000000000000..05d8123f3879341995e8898070bb6e1d22d91d2c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_concentration.md
@@ -0,0 +1 @@
+Process variable name of temperature.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_pressure.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..cadbd116efc42d12ab9e2bb570594c258aba6f3a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/process_variables/t_pressure.md
@@ -0,0 +1 @@
+Process variable name of pressure.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_decay_rate.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_decay_rate.md
new file mode 100644
index 0000000000000000000000000000000000000000..66473d364e13a2cbd454db05f258f25be4197899
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_decay_rate.md
@@ -0,0 +1 @@
+The decay rate \f$\vartheta\f$ used in the transport process.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid.md
new file mode 100644
index 0000000000000000000000000000000000000000..00aef4b95e01a85c9fb21ee4e2f36860cefd7dc1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid.md
@@ -0,0 +1 @@
+Defines fluid properties.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid_reference_density.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid_reference_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..7a9d865b220af527c3d672b644f8875d8b924c75
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_fluid_reference_density.md
@@ -0,0 +1 @@
+Parameter for the specification of the reference density of the fluid.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_molecular_diffusion_coefficient.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_molecular_diffusion_coefficient.md
new file mode 100644
index 0000000000000000000000000000000000000000..ba1c95fe91099233c51072710d9817f7dfeac8e8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_molecular_diffusion_coefficient.md
@@ -0,0 +1,2 @@
+Molecular diffusion coefficient \f$D_d\f$ used to compute the hydrodynamic
+dispersion tensor.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_retardation_factor.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_retardation_factor.md
new file mode 100644
index 0000000000000000000000000000000000000000..418ad86524af0595f76fcaa2a209f440e5227b33
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_retardation_factor.md
@@ -0,0 +1 @@
+Parameter for the specification of the retardation factor.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_longitudinal.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_longitudinal.md
new file mode 100644
index 0000000000000000000000000000000000000000..31b59b61b5f888522149c2bb4c47b192d41b84df
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_longitudinal.md
@@ -0,0 +1,2 @@
+The longitudinal dispersivity of chemical species \f$\beta_L\f$ used to compute
+the hydrodynamic dispersion tensor.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_transverse.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_transverse.md
new file mode 100644
index 0000000000000000000000000000000000000000..e4b2f6360e7d9053f5ce04455d4cc6ce022aab81
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_solute_dispersivity_transverse.md
@@ -0,0 +1,2 @@
+The transverse dispersivity of chemical species \f$\beta_T\f$ used to compute
+the hydrodynamic dispersion tensor.
diff --git a/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_specific_body_force.md
new file mode 100644
index 0000000000000000000000000000000000000000..36c8e1931c9001c35e937f9f2ac7014509f058ea
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RichardsComponentTransport/t_specific_body_force.md
@@ -0,0 +1,3 @@
+Specific body forces applied to fluid.
+
+It is usually used to apply gravitational forces.
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index aa8acef9f7287e25d49f6a699b533136a99b9d1b..75ec2adcc563fee3e91ac42a3a5796e6787a1216 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -4,6 +4,7 @@ APPEND_SOURCE_FILES(SOURCES)
 APPEND_SOURCE_FILES(SOURCES BoundaryCondition)
 APPEND_SOURCE_FILES(SOURCES CalculateSurfaceFlux)
 APPEND_SOURCE_FILES(SOURCES ComponentTransport)
+APPEND_SOURCE_FILES(SOURCES RichardsComponentTransport)
 APPEND_SOURCE_FILES(SOURCES Deformation)
 APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
 APPEND_SOURCE_FILES(SOURCES HeatConduction)
@@ -59,6 +60,7 @@ include(LIE/SmallDeformation/Tests.cmake)
 include(LiquidFlow/Tests.cmake)
 include(PhaseField/Tests.cmake)
 include(RichardsFlow/Tests.cmake)
+include(RichardsComponentTransport/Tests.cmake)
 include(SmallDeformation/Tests.cmake)
 include(TES/Tests.cmake)
 include(ThermoMechanics/Tests.cmake)
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index cbb55c4a4952bef43ce09d13201123fc5165fd2e..a237c8d51a025283d5ec63789de44132c0a1e035 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -306,6 +306,8 @@ public:
             auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
 
+            pos.setIntegrationPoint(ip);
+
             double T_int_pt = 0.0;
             double p_int_pt = 0.0;
             NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
diff --git a/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..99fe2891c179f6a579f5bdb3cf06dc7cc2a140de
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp
@@ -0,0 +1,113 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "CreatePorousMediaProperties.h"
+
+#include "BaseLib/reorderVector.h"
+
+#include "MaterialLib/PorousMedium/Permeability/createPermeabilityModel.h"
+#include "MaterialLib/PorousMedium/Porosity/createPorosityModel.h"
+#include "MaterialLib/PorousMedium/Storage/createStorageModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+
+#include "MeshLib/Mesh.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+PorousMediaProperties createPorousMediaProperties(
+    MeshLib::Mesh& mesh, BaseLib::ConfigTree const& porous_medium_configs)
+{
+    DBUG("Create PorousMediaProperties.");
+
+    std::vector<Eigen::MatrixXd> intrinsic_permeability_models;
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+        porosity_models;
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+        storage_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+        capillary_pressure_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+        relative_permeability_models;
+
+    std::vector<int> mat_ids;
+    for (auto const& porous_medium_config :
+         //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium}
+         porous_medium_configs.getConfigSubtreeList("porous_medium"))
+    {
+         //! \ogs_file_attr{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__id}
+        auto const id = porous_medium_config.getConfigAttribute<int>("id");
+        mat_ids.push_back(id);
+
+        auto const& porosity_config =
+            //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__porosity}
+            porous_medium_config.getConfigSubtree("porosity");
+        porosity_models.emplace_back(
+            MaterialLib::PorousMedium::createPorosityModel(porosity_config));
+
+        // Configuration for the intrinsic permeability
+        auto const& permeability_config =
+            //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__permeability}
+            porous_medium_config.getConfigSubtree("permeability");
+        intrinsic_permeability_models.emplace_back(
+            MaterialLib::PorousMedium::createPermeabilityModel(
+                permeability_config));
+
+        // Configuration for the specific storage.
+        auto const& storage_config =
+            //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__storage}
+            porous_medium_config.getConfigSubtree("storage");
+        storage_models.emplace_back(
+            MaterialLib::PorousMedium::createStorageModel(storage_config));
+
+        auto const& capillary_pressure_config =
+            //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__capillary_pressure}
+            porous_medium_config.getConfigSubtree("capillary_pressure");
+        auto capillary_pressure = MaterialLib::PorousMedium::createCapillaryPressureModel(
+            capillary_pressure_config);
+        capillary_pressure_models.emplace_back(std::move(capillary_pressure));
+
+        auto const& krel_config =
+            //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__relative_permeability}
+            porous_medium_config.getConfigSubtree("relative_permeability");
+        auto krel = MaterialLib::PorousMedium::createRelativePermeabilityModel(
+            krel_config);
+        relative_permeability_models.emplace_back(std::move(krel));
+    }
+
+    BaseLib::reorderVector(intrinsic_permeability_models, mat_ids);
+    BaseLib::reorderVector(porosity_models, mat_ids);
+    BaseLib::reorderVector(storage_models, mat_ids);
+
+    std::vector<int> material_ids(mesh.getNumberOfElements());
+    if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+    {
+        auto const& mesh_material_ids =
+            mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+        material_ids.reserve(mesh_material_ids->size());
+        std::copy(mesh_material_ids->cbegin(), mesh_material_ids->cend(),
+                  material_ids.begin());
+    }
+
+    return PorousMediaProperties{std::move(porosity_models),
+                                 std::move(intrinsic_permeability_models),
+                                 std::move(storage_models),
+                                 std::move(capillary_pressure_models),
+                                 std::move(relative_permeability_models),
+                                 std::move(material_ids)};
+}
+
+}  // namespace ComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..73e0ece1062340d303ccfb5839b813c8ebfe9e07
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "PorousMediaProperties.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+PorousMediaProperties createPorousMediaProperties(
+    MeshLib::Mesh& mesh, BaseLib::ConfigTree const& porous_media_config);
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..50a3377cd626089d55f2baebe075339a42309a3f
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
@@ -0,0 +1,153 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "CreateRichardsComponentTransportProcess.h"
+
+#include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
+
+#include "ProcessLib/Parameter/ConstantParameter.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "CreatePorousMediaProperties.h"
+#include "RichardsComponentTransportProcess.h"
+#include "RichardsComponentTransportProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+std::unique_ptr<Process> createRichardsComponentTransportProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__processes__process__type}
+    config.checkConfigParameter("type", "RichardsComponentTransport");
+
+    DBUG("Create RichardsComponentTransportProcess.");
+
+    // Process variable.
+
+    //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__process_variables}
+    auto const pv_config = config.getConfigSubtree("process_variables");
+
+    auto process_variables = findProcessVariables(
+        variables, pv_config,
+        {
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__process_variables__concentration}
+        "concentration",
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__process_variables__pressure}
+        "pressure"});
+
+    auto const& porous_medium_configs =
+        //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium}
+        config.getConfigSubtree("porous_medium");
+    PorousMediaProperties porous_media_properties{
+        createPorousMediaProperties(mesh, porous_medium_configs)};
+
+    //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__fluid}
+    auto const& fluid_config = config.getConfigSubtree("fluid");
+
+    auto fluid_properties =
+        MaterialLib::Fluid::createFluidProperties(fluid_config);
+
+    // Parameter for the density of the fluid.
+    auto& fluid_reference_density= findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__fluid_reference_density}
+        "fluid_reference_density", parameters, 1);
+    DBUG("Use \'%s\' as fluid_reference_density parameter.",
+         fluid_reference_density.name.c_str());
+
+    // Parameter for the longitudinal solute dispersivity.
+    auto const& molecular_diffusion_coefficient = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__molecular_diffusion_coefficient}
+        "molecular_diffusion_coefficient", parameters, 1);
+    DBUG("Use \'%s\' as molecular diffusion coefficient parameter.",
+         molecular_diffusion_coefficient.name.c_str());
+
+    // Parameter for the longitudinal solute dispersivity.
+    auto const& solute_dispersivity_longitudinal = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__solute_dispersivity_longitudinal}
+        "solute_dispersivity_longitudinal", parameters, 1);
+    DBUG("Use \'%s\' as longitudinal solute dispersivity parameter.",
+         solute_dispersivity_longitudinal.name.c_str());
+
+    // Parameter for the transverse solute dispersivity.
+    auto const& solute_dispersivity_transverse = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__solute_dispersivity_transverse}
+        "solute_dispersivity_transverse", parameters, 1);
+    DBUG("Use \'%s\' as transverse solute dispersivity parameter.",
+         solute_dispersivity_transverse.name.c_str());
+
+    // Parameter for the retardation factor.
+    auto const& retardation_factor = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__retardation_factor}
+        "retardation_factor", parameters, 1);
+
+    // Parameter for the decay rate.
+    auto const& decay_rate = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__decay_rate}
+        "decay_rate", parameters, 1);
+
+    // Specific body force parameter.
+    Eigen::VectorXd specific_body_force;
+    std::vector<double> const b =
+        //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__specific_body_force}
+        config.getConfigParameter<std::vector<double>>("specific_body_force");
+    assert(b.size() > 0 && b.size() < 4);
+    if (b.size() < mesh.getDimension())
+        OGS_FATAL(
+            "specific body force (gravity vector) has %d components, mesh "
+            "dimension is %d",
+            b.size(), mesh.getDimension());
+    bool const has_gravity = MathLib::toVector(b).norm() > 0;
+    if (has_gravity)
+    {
+        specific_body_force.resize(b.size());
+        std::copy_n(b.data(), b.size(), specific_body_force.data());
+    }
+
+    RichardsComponentTransportProcessData process_data{
+        std::move(porous_media_properties),
+        fluid_reference_density,
+        std::move(fluid_properties),
+        molecular_diffusion_coefficient,
+        solute_dispersivity_longitudinal,
+        solute_dispersivity_transverse,
+        retardation_factor,
+        decay_rate,
+        specific_body_force,
+        has_gravity};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"RichardsComponentTransport_concentration_pressure"});
+
+    ProcessLib::parseSecondaryVariables(config, secondary_variables,
+                                        named_function_caller);
+
+    return std::make_unique<RichardsComponentTransportProcess>(
+        mesh, std::move(jacobian_assembler), parameters, integration_order,
+        std::move(process_variables), std::move(process_data),
+        std::move(secondary_variables), std::move(named_function_caller));
+}
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..539b88eea044cd36b74b1d419e71e0d32b320bb7
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+std::unique_ptr<Process> createRichardsComponentTransportProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c89ce7c999ee640201bff1271a45bfe28a9c8679
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp
@@ -0,0 +1,57 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "PorousMediaProperties.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+int PorousMediaProperties::getMaterialID(SpatialPosition const& pos) const
+{
+    int const element_id = pos.getElementID().get();
+    return _material_ids[element_id];
+}
+
+MaterialLib::PorousMedium::Porosity const& PorousMediaProperties::getPorosity(
+    double /*t*/, SpatialPosition const& pos) const
+{
+    return *_porosity_models[getMaterialID(pos)];
+}
+
+Eigen::MatrixXd const& PorousMediaProperties::getIntrinsicPermeability(
+    double /*t*/, SpatialPosition const& pos) const
+{
+    return _intrinsic_permeability_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::Storage const&
+PorousMediaProperties::getSpecificStorage(double /*t*/,
+                                          SpatialPosition const& pos) const
+{
+    return *_specific_storage_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::CapillaryPressureSaturation const&
+PorousMediaProperties::getCapillaryPressureSaturationModel(
+    double /*t*/, SpatialPosition const& pos) const
+{
+    return *_capillary_pressure_saturation_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::RelativePermeability const&
+PorousMediaProperties::getRelativePermeability(
+    double /*t*/, SpatialPosition const& pos) const
+{
+    return *_relative_permeability_models[getMaterialID(pos)];
+}
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..6b538c77133f94dbb044f8cd60704ac86d9859fc
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h
@@ -0,0 +1,105 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+#include <Eigen/Dense>
+#include "BaseLib/reorderVector.h"
+
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
+
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+
+class PorousMediaProperties
+{
+public:
+    PorousMediaProperties(
+        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
+            porosity_models,
+        std::vector<Eigen::MatrixXd>&& intrinsic_permeability_models,
+        std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
+            specific_storage_models,
+        std::vector<std::unique_ptr<
+            MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+            capillary_pressure_saturation_models,
+        std::vector<
+            std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+            relative_permeability_models,
+        std::vector<int>&& material_ids)
+        : _porosity_models(std::move(porosity_models)),
+          _intrinsic_permeability_models(
+              std::move(intrinsic_permeability_models)),
+          _specific_storage_models(std::move(specific_storage_models)),
+          _capillary_pressure_saturation_models(
+              std::move(capillary_pressure_saturation_models)),
+          _relative_permeability_models(
+              std::move(relative_permeability_models)),
+          _material_ids(std::move(material_ids))
+    {
+    }
+
+    PorousMediaProperties(PorousMediaProperties&& other)
+        : _porosity_models(std::move(other._porosity_models)),
+          _intrinsic_permeability_models(other._intrinsic_permeability_models),
+          _specific_storage_models(std::move(other._specific_storage_models)),
+          _capillary_pressure_saturation_models(
+              std::move(other._capillary_pressure_saturation_models)),
+          _relative_permeability_models(
+              std::move(other._relative_permeability_models)),
+          _material_ids(other._material_ids)
+    {
+    }
+
+    MaterialLib::PorousMedium::Porosity const& getPorosity(
+        double t, SpatialPosition const& pos) const;
+
+    Eigen::MatrixXd const& getIntrinsicPermeability(
+        double t, SpatialPosition const& pos) const;
+
+    MaterialLib::PorousMedium::Storage const& getSpecificStorage(
+        double t, SpatialPosition const& pos) const;
+
+    MaterialLib::PorousMedium::CapillaryPressureSaturation const&
+    getCapillaryPressureSaturationModel(double t,
+                                        SpatialPosition const& pos) const;
+
+    MaterialLib::PorousMedium::RelativePermeability const&
+    getRelativePermeability(double t, SpatialPosition const& pos) const;
+
+private:
+    int getMaterialID(SpatialPosition const& pos) const;
+private:
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+        _porosity_models;
+    std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+        _specific_storage_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+        _capillary_pressure_saturation_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+        _relative_permeability_models;
+    std::vector<int> _material_ids;
+};
+
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..88c6c3f7bfe0256f4f70b4b06ed11e95e0026f06
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h
@@ -0,0 +1,344 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "RichardsComponentTransportFEM.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
+    LocalAssemblerData(
+        MeshLib::Element const& element,
+        std::size_t const local_matrix_size,
+        bool is_axially_symmetric,
+        unsigned const integration_order,
+        RichardsComponentTransportProcessData const& process_data)
+    : _element_id(element.getID()),
+      _process_data(process_data),
+      _integration_method(integration_order)
+{
+    // This assertion is valid only if all nodal d.o.f. use the same shape
+    // matrices.
+    assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+    (void)local_matrix_size;
+
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+    _ip_data.reserve(n_integration_points);
+
+    auto const shape_matrices =
+        initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
+                          GlobalDim>(element, is_axially_symmetric,
+                                     _integration_method);
+
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
+    {
+        auto const& sm = shape_matrices[ip];
+        double const integration_factor = sm.integralMeasure * sm.detJ;
+        _ip_data.emplace_back(
+            sm.N, sm.dNdx,
+            _integration_method.getWeightedPoint(ip).getWeight() *
+                integration_factor,
+            sm.N.transpose() * sm.N * integration_factor *
+                _integration_method.getWeightedPoint(ip).getWeight());
+    }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::assemble(
+    double const t,
+    std::vector<double> const& local_x,
+    std::vector<double>& local_M_data,
+    std::vector<double>& local_K_data,
+    std::vector<double>& local_b_data)
+{
+    auto const local_matrix_size = local_x.size();
+    // This assertion is valid only if all nodal d.o.f. use the same shape
+    // matrices.
+    assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+    auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+        local_M_data, local_matrix_size, local_matrix_size);
+    auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+        local_K_data, local_matrix_size, local_matrix_size);
+    auto local_b = MathLib::createZeroedVector<LocalVectorType>(
+        local_b_data, local_matrix_size);
+
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    SpatialPosition pos;
+    pos.setElementID(_element_id);
+
+    auto p_nodal_values = Eigen::Map<const NodalVectorType>(
+        &local_x[pressure_index], pressure_size);
+
+    auto const& b = _process_data.specific_body_force;
+
+    MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+    GlobalDimMatrixType const& I(
+        GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
+
+    auto KCC = local_K.template block<concentration_size, concentration_size>(
+        concentration_index, concentration_index);
+    auto MCC = local_M.template block<concentration_size, concentration_size>(
+        concentration_index, concentration_index);
+    auto Kpp = local_K.template block<pressure_size, pressure_size>(
+        pressure_index, pressure_index);
+    auto Mpp = local_M.template block<pressure_size, pressure_size>(
+        pressure_index, pressure_index);
+    auto Bp = local_b.template block<pressure_size, 1>(pressure_index, 0);
+
+    for (std::size_t ip(0); ip < n_integration_points; ++ip)
+    {
+        pos.setIntegrationPoint(ip);
+
+        // \todo the argument to getValue() has to be changed for non
+        // constant storage model
+        auto const specific_storage =
+            _process_data.porous_media_properties.getSpecificStorage(t, pos)
+                .getValue(0.0);
+
+        auto const& ip_data = _ip_data[ip];
+        auto const& N = ip_data.N;
+        auto const& dNdx = ip_data.dNdx;
+        auto const& w = ip_data.integration_weight;
+
+        double C_int_pt = 0.0;
+        double p_int_pt = 0.0;
+        // Order matters: First C, then p!
+        NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
+
+        double const pc_int_pt = -p_int_pt;
+        double const Sw = _process_data.porous_media_properties
+                              .getCapillaryPressureSaturationModel(t, pos)
+                              .getSaturation(pc_int_pt);
+
+        double const dSw_dpc =
+            1. / _process_data.porous_media_properties
+                     .getCapillaryPressureSaturationModel(t, pos)
+                     .getdPcdS(Sw);
+
+        // \todo the first argument has to be changed for non constant
+        // porosity model
+        auto const porosity =
+            _process_data.porous_media_properties.getPorosity(t, pos).getValue(
+                0.0, C_int_pt);
+
+        auto const retardation_factor =
+            _process_data.retardation_factor(t, pos)[0];
+
+        auto const& solute_dispersivity_transverse =
+            _process_data.solute_dispersivity_transverse(t, pos)[0];
+        auto const& solute_dispersivity_longitudinal =
+            _process_data.solute_dispersivity_longitudinal(t, pos)[0];
+
+        // Use the fluid density model to compute the density
+        vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::C)] =
+            C_int_pt;
+        vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
+            p_int_pt;
+        auto const density = _process_data.fluid_properties->getValue(
+            MaterialLib::Fluid::FluidPropertyType::Density, vars);
+        auto const& decay_rate = _process_data.decay_rate(t, pos)[0];
+        auto const& molecular_diffusion_coefficient =
+            _process_data.molecular_diffusion_coefficient(t, pos)[0];
+
+        auto const& K =
+            _process_data.porous_media_properties.getIntrinsicPermeability(t,
+                                                                           pos);
+        auto const& k_rel = _process_data.porous_media_properties
+                                .getRelativePermeability(t, pos)
+                                .getValue(Sw);
+        // Use the viscosity model to compute the viscosity
+        auto const mu = _process_data.fluid_properties->getValue(
+            MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+        auto const K_times_k_rel_over_mu = K * (k_rel / mu);
+
+        GlobalDimVectorType const velocity =
+            _process_data.has_gravity
+                ? GlobalDimVectorType(-K_times_k_rel_over_mu *
+                                      (dNdx * p_nodal_values - density * b))
+                : GlobalDimVectorType(-K_times_k_rel_over_mu * dNdx *
+                                      p_nodal_values);
+
+        double const velocity_magnitude = velocity.norm();
+        GlobalDimMatrixType const hydrodynamic_dispersion =
+            velocity_magnitude != 0.0
+                ? GlobalDimMatrixType(
+                      (porosity * molecular_diffusion_coefficient +
+                       solute_dispersivity_transverse * velocity_magnitude) *
+                          I +
+                      (solute_dispersivity_longitudinal -
+                       solute_dispersivity_transverse) /
+                          velocity_magnitude * velocity * velocity.transpose())
+                : GlobalDimMatrixType(
+                      (porosity * molecular_diffusion_coefficient +
+                       solute_dispersivity_transverse * velocity_magnitude) *
+                      I);
+
+        // matrix assembly
+        KCC.noalias() +=
+            (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
+             N.transpose() * velocity.transpose() * dNdx +
+             N.transpose() * decay_rate * porosity * retardation_factor * N) *
+            w;
+        MCC.noalias() += w * N.transpose() * porosity * retardation_factor * N;
+        Kpp.noalias() += w * dNdx.transpose() * K_times_k_rel_over_mu * dNdx;
+        // \TODO Extend to pressure dependent density.
+        double const drhow_dp(0.0);
+        Mpp.noalias() += (specific_storage * Sw + porosity * Sw * drhow_dp -
+                          porosity * dSw_dpc) *
+                         ip_data.mass_operator;
+
+        if (_process_data.has_gravity)
+            Bp += w * density * dNdx.transpose() * K_times_k_rel_over_mu * b;
+        /* with Oberbeck-Boussing assumption density difference only exists
+         * in buoyancy effects */
+    }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+std::vector<double> const&
+LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
+    getIntPtDarcyVelocity(const double t,
+                          GlobalVector const& current_solution,
+                          NumLib::LocalToGlobalIndexMap const& dof_table,
+                          std::vector<double>& cache) const
+{
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    auto const indices = NumLib::getIndices(_element_id, dof_table);
+    assert(!indices.empty());
+    auto const local_x = current_solution.get(indices);
+
+    cache.clear();
+    auto cache_mat = MathLib::createZeroedMatrix<
+        Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+        cache, GlobalDim, n_integration_points);
+
+    SpatialPosition pos;
+    pos.setElementID(_element_id);
+
+    MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+    auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+        &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+
+    for (unsigned ip = 0; ip < n_integration_points; ++ip)
+    {
+        auto const& ip_data = _ip_data[ip];
+        auto const& N = ip_data.N;
+        auto const& dNdx = ip_data.dNdx;
+
+        pos.setIntegrationPoint(ip);
+
+        auto const& K =
+            _process_data.porous_media_properties.getIntrinsicPermeability(t,
+                                                                           pos);
+        auto const mu = _process_data.fluid_properties->getValue(
+            MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+
+        double C_int_pt = 0.0;
+        double p_int_pt = 0.0;
+        NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
+
+        // saturation
+        double const pc_int_pt = -p_int_pt;
+        double const Sw = _process_data.porous_media_properties
+                              .getCapillaryPressureSaturationModel(t, pos)
+                              .getSaturation(pc_int_pt);
+
+        auto const& k_rel = _process_data.porous_media_properties
+                                .getRelativePermeability(t, pos)
+                                .getValue(Sw);
+
+        cache_mat.col(ip).noalias() = -dNdx * p_nodal_values;
+        if (_process_data.has_gravity)
+        {
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+
+            auto const rho_w = _process_data.fluid_properties->getValue(
+                MaterialLib::Fluid::FluidPropertyType::Density, vars);
+            auto const b = _process_data.specific_body_force;
+            // here it is assumed that the vector b is directed 'downwards'
+            cache_mat.col(ip).noalias() += rho_w * b;
+        }
+        cache_mat.col(ip).noalias() = k_rel / mu * (K * cache_mat.col(ip));
+    }
+    return cache;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+Eigen::Map<const Eigen::RowVectorXd>
+LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::getShapeMatrix(
+    const unsigned integration_point) const
+{
+    auto const& N = _ip_data[integration_point].N;
+
+    // assumes N is stored contiguously in memory
+    return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+std::vector<double> const&
+LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
+    getIntPtSaturation(const double t,
+                       GlobalVector const& current_solution,
+                       NumLib::LocalToGlobalIndexMap const& dof_table,
+                       std::vector<double>& cache) const
+{
+    SpatialPosition pos;
+    pos.setElementID(_element_id);
+
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    auto const indices = NumLib::getIndices(_element_id, dof_table);
+    assert(!indices.empty());
+    auto const local_x = current_solution.get(indices);
+
+    cache.clear();
+    auto cache_vec = MathLib::createZeroedVector<LocalVectorType>(
+        cache, n_integration_points);
+
+    for (unsigned ip = 0; ip < n_integration_points; ++ip)
+    {
+        auto const& ip_data = _ip_data[ip];
+        auto const& N = ip_data.N;
+
+        double C_int_pt = 0.0;
+        double p_int_pt = 0.0;
+        NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
+
+        // saturation
+        double const pc_int_pt = -p_int_pt;
+        double const Sw = _process_data.porous_media_properties
+                              .getCapillaryPressureSaturationModel(t, pos)
+                              .getSaturation(pc_int_pt);
+        cache_vec[ip] = Sw;
+    }
+
+    return cache;
+}
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c51bd8e93e63454a808ca7eb4e91d6af5490574
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
@@ -0,0 +1,146 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+#include "RichardsComponentTransportProcessData.h"
+#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType,
+          typename NodalMatrixType>
+struct IntegrationPointData final
+{
+    IntegrationPointData(NodalRowVectorType const& N_,
+                         GlobalDimNodalMatrixType const& dNdx_,
+                         double const& integration_weight_,
+                         NodalMatrixType const mass_operator_)
+        : N(N_),
+          dNdx(dNdx_),
+          integration_weight(integration_weight_),
+          mass_operator(mass_operator_)
+    {
+    }
+
+    NodalRowVectorType const N;
+    GlobalDimNodalMatrixType const dNdx;
+    double const integration_weight;
+    NodalMatrixType const mass_operator;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+const unsigned NUM_NODAL_DOF = 2;
+
+class RichardsComponentTransportLocalAssemblerInterface
+    : public ProcessLib::LocalAssemblerInterface,
+      public NumLib::ExtrapolatableElement
+{
+public:
+    virtual std::vector<double> const& getIntPtSaturation(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class LocalAssemblerData
+    : public RichardsComponentTransportLocalAssemblerInterface
+{
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+    using LocalMatrixType = typename ShapeMatricesType::template MatrixType<
+        NUM_NODAL_DOF * ShapeFunction::NPOINTS,
+        NUM_NODAL_DOF * ShapeFunction::NPOINTS>;
+    using LocalVectorType =
+        typename ShapeMatricesType::template VectorType<NUM_NODAL_DOF *
+                                                        ShapeFunction::NPOINTS>;
+
+    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+    using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+
+    using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+    using GlobalDimNodalMatrixType =
+        typename ShapeMatricesType::GlobalDimNodalMatrixType;
+    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+    using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+
+public:
+    LocalAssemblerData(
+        MeshLib::Element const& element,
+        std::size_t const local_matrix_size,
+        bool is_axially_symmetric,
+        unsigned const integration_order,
+        RichardsComponentTransportProcessData const& process_data);
+
+    void assemble(double const t, std::vector<double> const& local_x,
+                  std::vector<double>& local_M_data,
+                  std::vector<double>& local_K_data,
+                  std::vector<double>& local_b_data) override;
+
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override;
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override;
+
+    std::vector<double> const& getIntPtSaturation(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override;
+private:
+    unsigned const _element_id;
+    RichardsComponentTransportProcessData const& _process_data;
+
+    IntegrationMethod const _integration_method;
+    std::vector<
+        IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType,
+                             NodalMatrixType>,
+        Eigen::aligned_allocator<IntegrationPointData<
+            NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
+        _ip_data;
+
+    static const int concentration_index = 0;
+    static const int concentration_size = ShapeFunction::NPOINTS;
+    static const int pressure_index = ShapeFunction::NPOINTS;
+    static const int pressure_size = ShapeFunction::NPOINTS;
+};
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
+
+#include "RichardsComponentTransportFEM-impl.h"
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..67a93c25251d40ca8553e8b6f500677980f454c6
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
@@ -0,0 +1,92 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "RichardsComponentTransportProcess.h"
+
+#include <cassert>
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+RichardsComponentTransportProcess::RichardsComponentTransportProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    RichardsComponentTransportProcessData&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller)
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              integration_order, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller)),
+      _process_data(std::move(process_data))
+{
+}
+
+void RichardsComponentTransportProcess::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
+    ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+        mesh.getDimension(), mesh.getElements(), dof_table,
+        pv.getShapeFunctionOrder(), _local_assemblers,
+        mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+    _secondary_variables.addSecondaryVariable(
+        "darcy_velocity",
+        makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+                         _local_assemblers,
+                         &RichardsComponentTransportLocalAssemblerInterface::
+                             getIntPtDarcyVelocity));
+
+    _secondary_variables.addSecondaryVariable(
+        "saturation",
+        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+                         &RichardsComponentTransportLocalAssemblerInterface::
+                             getIntPtSaturation));
+}
+
+void RichardsComponentTransportProcess::assembleConcreteProcess(
+    const double t,
+    GlobalVector const& x,
+    GlobalMatrix& M,
+    GlobalMatrix& K,
+    GlobalVector& b)
+{
+    DBUG("Assemble RichardsComponentTransportProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+        *_local_to_global_index_map, t, x, M, K, b, _coupling_term);
+}
+
+void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
+    const double t, GlobalVector const& x, GlobalVector const& xdot,
+    const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b, GlobalMatrix& Jac)
+{
+    DBUG("AssembleWithJacobian RichardsComponentTransportProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+        dx_dx, M, K, b, Jac, _coupling_term);
+}
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
+
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..9e81aac5a3bbd297314d6e95e0fd3a22521cb878
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
@@ -0,0 +1,145 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "RichardsComponentTransportFEM.h"
+#include "RichardsComponentTransportProcessData.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+/**
+ * # RichardsComponentTransport process
+ *
+ * ## Governing Equations
+ *
+ * ### Richards Flow
+ *
+ * The flow process is described by
+ * \f[
+ * \phi \frac{\partial \rho_w}{\partial p} \frac{\partial p}{\partial t} S
+ *     - \phi \rho_w \frac{\partial S}{\partial p_c}
+ *       \frac{\partial p_c}{\partial t}
+ *     - \nabla \cdot \left[\rho_w \frac{k_{\mathrm{rel}} \kappa}{\mu}
+ *       \nabla \left( p + \rho_w g z \right)\right]
+ *     - Q_p = 0,
+ * \f]
+ * where
+ * - \f$\phi\f$ is the porosity,
+ * - \f$S\f$ is the saturation,
+ * - \f$p\f$ is the pressure,
+ * - \f$k_{\mathrm{rel}}\f$ is the relativ permeability (depending on \f$S\f$),
+ * - \f$\kappa\f$ is the specific permeability,
+ * - \f$\mu\f$ is viscosity of the fluid,
+ * - \f$\rho_w\f$ is the mass density of the fluid, and
+ * - \f$g\f$ is the gravitational acceleration.
+ *
+ * Here it is assumed, that
+ * - the porosity is constant and
+ * - the pressure of the gas phase is zero.
+ *
+ * The capillary pressure is given by
+ * \f[
+ * p_c = \frac{\rho_w g}{\alpha}
+ *  \left[S_{\mathrm{eff}}^{-\frac{1}{m}} - 1\right]^{\frac{1}{n}}
+ * \f]
+ * and the effective saturation by
+ * \f[
+ * S_{\mathrm{eff}} = \frac{S-S_r}{S_{\max} - S_r}
+ * \f]
+ *
+ * ### Mass Transport
+ * The mass transport process is described by
+ * \f[
+ * \phi R \frac{\partial C}{\partial t}
+    + \nabla \cdot \left(\vec{q} C - D \nabla C \right)
+        + \phi R \vartheta C - Q_C = 0
+ * \f]
+ * where
+ * - \f$R\f$ is the retardation factor,
+ * - \f$C\f$ is the concentration,
+ * - \f$\vec{q} = \frac{k_{\mathrm{rel}} \kappa}{\mu(C)}
+ *   \nabla \left( p + \rho_w g z \right)\f$ is the Darcy velocity,
+ * - \f$D\f$ is the hydrodynamic dispersion tensor,
+ * - \f$\vartheta\f$ is the decay rate.
+ *
+ * For the hydrodynamic dispersion tensor the relation
+ * \f[
+ * D = (\phi D_d + \beta_T \|\vec{q}\|) I + (\beta_L - \beta_T) \frac{\vec{q}
+ * \vec{q}^T}{\|\vec{q}\|}
+ * \f]
+ * is implemented, where \f$D_d\f$ is the molecular diffusion coefficient,
+ * \f$\beta_L\f$ the longitudinal dispersivity of chemical species, and
+ * \f$\beta_T\f$ the transverse dispersivity of chemical species.
+ *
+ * The implementation uses a monolithic approach, i.e., both processes
+ * are assembled within one global system of equations.
+ *
+ * ## Process Coupling
+ *
+ * The advective term of the concentration equation is given by the Richards
+ * flow process, i.e., the concentration distribution depends on
+ * darcy velocity of the Richards flow process. On the other hand the
+ * concentration dependencies of the viscosity and density in the groundwater
+ * flow couples the unsaturated H process to the C process.
+ *
+ * \note At the moment there is not any coupling by source or sink terms, i.e.,
+ * the coupling is implemented only by density changes due to concentration
+ * changes in the buoyance term in the groundwater flow. This coupling schema is
+ * referred to as the Boussinesq approximation.
+ * */
+class RichardsComponentTransportProcess final : public Process
+{
+public:
+    RichardsComponentTransportProcess(
+        MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        unsigned const integration_order,
+        std::vector<std::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        RichardsComponentTransportProcessData&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller);
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override { return false; }
+    //! @}
+
+private:
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
+
+    void assembleWithJacobianConcreteProcess(
+        const double t, GlobalVector const& x, GlobalVector const& xdot,
+        const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+    RichardsComponentTransportProcessData _process_data;
+
+    std::vector<
+        std::unique_ptr<RichardsComponentTransportLocalAssemblerInterface>>
+        _local_assemblers;
+};
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..97c4c58a803b9714ecd38469f84cdf201a4036c8
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h
@@ -0,0 +1,94 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+#include "PorousMediaProperties.h"
+
+namespace ProcessLib
+{
+template <typename ReturnType>
+struct Parameter;
+
+namespace RichardsComponentTransport
+{
+struct RichardsComponentTransportProcessData
+{
+    RichardsComponentTransportProcessData(
+        PorousMediaProperties&& porous_media_properties_,
+        ProcessLib::Parameter<double> const& fluid_reference_density_,
+        std::unique_ptr<MaterialLib::Fluid::FluidProperties>&&
+            fluid_properties_,
+        ProcessLib::Parameter<double> const& molecular_diffusion_coefficient_,
+        ProcessLib::Parameter<double> const& solute_dispersivity_longitudinal_,
+        ProcessLib::Parameter<double> const& solute_dispersivity_transverse_,
+        ProcessLib::Parameter<double> const& retardation_factor_,
+        ProcessLib::Parameter<double> const& decay_rate_,
+        Eigen::VectorXd const& specific_body_force_,
+        bool const has_gravity_)
+        : porous_media_properties(std::move(porous_media_properties_)),
+          fluid_reference_density(fluid_reference_density_),
+          fluid_properties(std::move(fluid_properties_)),
+          molecular_diffusion_coefficient(molecular_diffusion_coefficient_),
+          solute_dispersivity_longitudinal(solute_dispersivity_longitudinal_),
+          solute_dispersivity_transverse(solute_dispersivity_transverse_),
+          retardation_factor(retardation_factor_),
+          decay_rate(decay_rate_),
+          specific_body_force(specific_body_force_),
+          has_gravity(has_gravity_)
+    {
+    }
+
+    RichardsComponentTransportProcessData(
+        RichardsComponentTransportProcessData&& other)
+        : porous_media_properties(std::move(other.porous_media_properties)),
+          fluid_reference_density(other.fluid_reference_density),
+          fluid_properties(other.fluid_properties.release()),
+          molecular_diffusion_coefficient(
+              other.molecular_diffusion_coefficient),
+          solute_dispersivity_longitudinal(
+              other.solute_dispersivity_longitudinal),
+          solute_dispersivity_transverse(other.solute_dispersivity_transverse),
+          retardation_factor(other.retardation_factor),
+          decay_rate(other.decay_rate),
+          specific_body_force(other.specific_body_force),
+          has_gravity(other.has_gravity)
+    {
+    }
+
+    //! Copies are forbidden.
+    RichardsComponentTransportProcessData(
+        RichardsComponentTransportProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(RichardsComponentTransportProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(RichardsComponentTransportProcessData&&) = delete;
+
+    PorousMediaProperties porous_media_properties;
+    Parameter<double> const& fluid_reference_density;
+    std::unique_ptr<MaterialLib::Fluid::FluidProperties> fluid_properties;
+    Parameter<double> const& molecular_diffusion_coefficient;
+    Parameter<double> const& solute_dispersivity_longitudinal;
+    Parameter<double> const& solute_dispersivity_transverse;
+    Parameter<double> const& retardation_factor;
+    Parameter<double> const& decay_rate;
+    Eigen::VectorXd const specific_body_force;
+    bool const has_gravity;
+};
+
+}  // namespace RichardsComponentTransport
+}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/Tests.cmake b/ProcessLib/RichardsComponentTransport/Tests.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..09139183c5da0aa1e6c2a8bfb4965340682956ba
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/Tests.cmake
@@ -0,0 +1,292 @@
+AddTest(
+    NAME RichardsComponentTransport_1D_Padilla_NaCl1
+    PATH Parabolic/RichardsComponentTransport/Padilla/
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS Padilla_NaCl1.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    ABSTOL 1e-2 RELTOL 1e-10 # darcy velocity has a large absolute error
+    DIFF_DATA
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu concentration concentration
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu pressure pressure
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl1_pcs_0_ts_0_t_0.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu Padilla_NaCl1_pcs_0_ts_10005_t_20000.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu Padilla_NaCl1_pcs_0_ts_1050_t_2090.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu Padilla_NaCl1_pcs_0_ts_2050_t_4090.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu Padilla_NaCl1_pcs_0_ts_2100_t_4190.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu Padilla_NaCl1_pcs_0_ts_2150_t_4290.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu Padilla_NaCl1_pcs_0_ts_2200_t_4390.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu Padilla_NaCl1_pcs_0_ts_2250_t_4490.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu Padilla_NaCl1_pcs_0_ts_2300_t_4590.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu Padilla_NaCl1_pcs_0_ts_2350_t_4690.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu Padilla_NaCl1_pcs_0_ts_2400_t_4790.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu Padilla_NaCl1_pcs_0_ts_2450_t_4890.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu Padilla_NaCl1_pcs_0_ts_2500_t_4990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu Padilla_NaCl1_pcs_0_ts_2550_t_5090.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu Padilla_NaCl1_pcs_0_ts_2600_t_5190.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu Padilla_NaCl1_pcs_0_ts_2650_t_5290.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu Padilla_NaCl1_pcs_0_ts_2700_t_5390.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu Padilla_NaCl1_pcs_0_ts_2750_t_5490.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu Padilla_NaCl1_pcs_0_ts_2800_t_5590.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu Padilla_NaCl1_pcs_0_ts_2850_t_5690.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu Padilla_NaCl1_pcs_0_ts_2900_t_5790.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu Padilla_NaCl1_pcs_0_ts_2950_t_5890.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl1_pcs_0_ts_3000_t_5990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu Padilla_NaCl1_pcs_0_ts_3050_t_6090.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu Padilla_NaCl1_pcs_0_ts_3100_t_6190.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu Padilla_NaCl1_pcs_0_ts_3150_t_6290.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu Padilla_NaCl1_pcs_0_ts_3200_t_6390.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu Padilla_NaCl1_pcs_0_ts_3500_t_6990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl1_pcs_0_ts_4500_t_8990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl1_pcs_0_ts_5500_t_10990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu Padilla_NaCl1_pcs_0_ts_550_t_1090.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl1_pcs_0_ts_6500_t_12990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl1_pcs_0_ts_7500_t_14990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl1_pcs_0_ts_8500_t_16990.000000.vtu saturation saturation
+    Padilla_NaCl1/Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl1_pcs_0_ts_9500_t_18990.000000.vtu saturation saturation
+)
+
+AddTest(
+    NAME RichardsComponentTransport_1D_Padilla_NaCl6
+    PATH Parabolic/RichardsComponentTransport/Padilla/
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS Padilla_NaCl6.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    ABSTOL 1e-2 RELTOL 1e-10 # darcy velocity has a large absolute error
+    DIFF_DATA
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu concentration concentration
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu darcy_velocity darcy_velocity
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu pressure pressure
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu Padilla_NaCl6_pcs_0_ts_0_t_0.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu Padilla_NaCl6_pcs_0_ts_10000_t_19990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu Padilla_NaCl6_pcs_0_ts_1000_t_1990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu Padilla_NaCl6_pcs_0_ts_10500_t_20990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu Padilla_NaCl6_pcs_0_ts_11000_t_21990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu Padilla_NaCl6_pcs_0_ts_11500_t_22990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu Padilla_NaCl6_pcs_0_ts_12000_t_23990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu Padilla_NaCl6_pcs_0_ts_12500_t_24990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu Padilla_NaCl6_pcs_0_ts_13000_t_25990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu Padilla_NaCl6_pcs_0_ts_13500_t_26990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu Padilla_NaCl6_pcs_0_ts_14000_t_27990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu Padilla_NaCl6_pcs_0_ts_15000_t_29990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu Padilla_NaCl6_pcs_0_ts_16000_t_31990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu Padilla_NaCl6_pcs_0_ts_17000_t_33990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu Padilla_NaCl6_pcs_0_ts_18000_t_35990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu Padilla_NaCl6_pcs_0_ts_19000_t_37990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu Padilla_NaCl6_pcs_0_ts_20000_t_39990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu Padilla_NaCl6_pcs_0_ts_20005_t_40000.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu Padilla_NaCl6_pcs_0_ts_2000_t_3990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu Padilla_NaCl6_pcs_0_ts_3000_t_5990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu Padilla_NaCl6_pcs_0_ts_4000_t_7990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu Padilla_NaCl6_pcs_0_ts_4500_t_8990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu Padilla_NaCl6_pcs_0_ts_5000_t_9990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu Padilla_NaCl6_pcs_0_ts_5500_t_10990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu Padilla_NaCl6_pcs_0_ts_6000_t_11990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu Padilla_NaCl6_pcs_0_ts_6500_t_12990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu Padilla_NaCl6_pcs_0_ts_7000_t_13990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu Padilla_NaCl6_pcs_0_ts_7500_t_14990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu Padilla_NaCl6_pcs_0_ts_8000_t_15990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu Padilla_NaCl6_pcs_0_ts_8500_t_16990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu Padilla_NaCl6_pcs_0_ts_9000_t_17990.000000.vtu saturation saturation
+    Padilla_NaCl6/Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu Padilla_NaCl6_pcs_0_ts_9500_t_18990.000000.vtu saturation saturation
+)
+
diff --git a/web/content/docs/benchmarks/richards-flow/Padilla_NaCl1_poreVolume-conc.csv b/web/content/docs/benchmarks/richards-flow/Padilla_NaCl1_poreVolume-conc.csv
new file mode 100644
index 0000000000000000000000000000000000000000..bd123086e7c0f2766f3a31f7199b0e74666e9796
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-flow/Padilla_NaCl1_poreVolume-conc.csv
@@ -0,0 +1,34 @@
+x, y
+0.85, -0.01
+0.87, -0.01
+0.89, 0
+0.91, 0.02
+0.92, 0.04
+0.93, 0.07
+0.94, 0.11
+0.95, 0.16
+0.96, 0.2
+0.97, 0.25
+0.97, 0.31
+0.98, 0.37
+0.99, 0.44
+0.99, 0.49
+1, 0.54
+1.01, 0.61
+1.02, 0.67
+1.03, 0.76
+1.04, 0.8
+1.05, 0.86
+1.07, 0.91
+1.09, 0.95
+1.1, 0.97
+1.13, 0.99
+1.15, 0.99
+1.19, 0.99
+1.21, 0.99
+1.23, 0.99
+1.25, 0.99
+1.29, 1
+1.36, 1
+1.4, 1
+1.43, 1
diff --git a/web/content/docs/benchmarks/richards-flow/Padilla_NaCl6_poreVolume-conc.csv b/web/content/docs/benchmarks/richards-flow/Padilla_NaCl6_poreVolume-conc.csv
new file mode 100644
index 0000000000000000000000000000000000000000..72344dae429bbf7da14d29f1cfcbd6262c0dfb01
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-flow/Padilla_NaCl6_poreVolume-conc.csv
@@ -0,0 +1,45 @@
+x, y
+0.64, -0.01
+0.67, -0.01
+0.7, -0.01
+0.73, 0.01
+0.76, 0.04
+0.77, 0.08
+0.8, 0.12
+0.82, 0.16
+0.84, 0.21
+0.87, 0.28
+0.87, 0.32
+0.9, 0.37
+0.92, 0.43
+0.94, 0.47
+0.97, 0.53
+0.98, 0.57
+1.01, 0.62
+1.05, 0.68
+1.07, 0.71
+1.1, 0.74
+1.12, 0.77
+1.14, 0.79
+1.16, 0.8
+1.18, 0.82
+1.2, 0.84
+1.24, 0.88
+1.26, 0.89
+1.28, 0.9
+1.3, 0.91
+1.32, 0.92
+1.35, 0.93
+1.39, 0.94
+1.43, 0.95
+1.45, 0.96
+1.49, 0.96
+1.52, 0.97
+1.56, 0.97
+1.6, 0.98
+1.64, 0.98
+1.69, 0.98
+1.72, 0.98
+1.75, 0.98
+1.82, 0.98
+1.85, 0.99
diff --git a/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Equations.pdf b/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Equations.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..2a5d8703d8207ca84d057b60385a8c42f129d753
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Equations.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bd5a9e1ffabcca2e9faf2e8d99e50a4c16c63646a8038acd5bd8a885730f63a0
+size 296677
diff --git a/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Padilla.png b/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Padilla.png
new file mode 100644
index 0000000000000000000000000000000000000000..dbb2ddad72192d251983b6bbc04cd79fe34d9f4d
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-flow/RichardsComponentTransport_Padilla.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:48c8aea7742c981d58c6d6afb9a0c35c322a927ad441f7011f9fc47a3008ab65
+size 31188
diff --git a/web/content/docs/benchmarks/richards-flow/richards-component-transport.md b/web/content/docs/benchmarks/richards-flow/richards-component-transport.md
new file mode 100644
index 0000000000000000000000000000000000000000..76d821f75d54afa3264800f994bf483b861d5a06
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-flow/richards-component-transport.md
@@ -0,0 +1,55 @@
++++
+author = "Marc Walther"
+weight = 142
+project = "Parabolic/RichardsComponentTransport/Padilla/unsaturated_NaCl6.prj"
+date = "2017-08-25T14:41:09+01:00"
+title = "Unsaturated Mass Transport"
+
+[menu]
+  [menu.benchmarks]
+    parent = "richards-flow"
+
++++
+
+{{< project-link >}}
+
+
+## Overview
+
+The Richards equation is often used to describe water movement in the unsaturated zone, while the mass transport equation describes solute movement in the liquid phase. Here, we use a monolithic approach to simulate mass transport in an unsaturated medium.
+
+The development of the equation system is given in [this PDF](../RichardsComponentTransport_Equations.pdf). In the following, we present a numerical benchmark that uses experimental data as reference.
+
+
+## Problem description
+
+We use the Padilla et al. (1999) problem which features a series of saturated and unsaturated laboratory column experiments to validate the implementation. In specific, we refer to experiments "NaCl1" (saturated flow) and "NaCl6" (unsaturated flow) - see also table 1 in Padilla et al. (1999) - and use the following parameters as model input.
+
+
+### Model setup
+
+We use a 1D domain with $0 < x < 0.25 m$, and a spatial resolution of $1.25 mm$ with a total of 200 elements. Top boundary conditions are concentration $c = 1$ (as Dirichlet) for both setups, and a specific flux $q\_{NaCl1} = 2.12789E-05$ and $q\_{NaCl6} = 6.46004E-06$ (as Neumann) for saturated and unsaturated cases, respectively. Bottom boundary conditions are a free exit boundary for mass transport, and Dirichlet pressure values are chosen so that the saturation in the column follows the respective values according to table 1 in Padella et al. (1999) as $p\_{NaCl1} = 0 Pa$ and $p\_{NaCl6} = -4800 Pa$.
+
+Saturated intrinsic permeability is calculated from given flux and pressure gradient information as $\kappa = 1.174E-10 m^2 $; total porosity is $\theta = 0.45$; van-Genuchten-values are $m = 0.789$ (from $n = 4.74$ via $m = 1-1/n$), bubbling pressure $p_d = 3633.33 Pa$ (from $\alpha = 0.027 cm^{-1}$ via $p_c = \rho \cdot g / \alpha$), residual saturation $s_r = 0.1689$. Molecular diffusion coefficient was set to $D_m = 1e-9 m^2/s$; dispersivities were chosen from table 2 in Padilla et al. (1999) and set to $\alpha\_{NaCl1} = 3.4173E-04 m$ and $\alpha\_{NaCl6} = 4.6916E-03 m$ for saturated and unsaturated conditions, respectively.
+
+Initial conditions are $c = 0$ and hydrostatic pressure conditions with steady state flow for both scenarios.
+
+[The NaCl1 project file is here. ](../../../../../Tests/Data/Parabolic/RichardsComponentTransport/Padilla_NaCl1.prj)
+[The NaCl6 project file is here. ](../../../../../Tests/Data/Parabolic/RichardsComponentTransport/Padilla_NaCl6.prj)
+
+
+### Results
+
+The figure below shows breakthrough curves vs experimental result at the bottom of the simulation domain, together with averaged saturation values at the two observation points with distance of 0.075 cm from both ends of the column (as stated in Padilla et al., 1999) over pore volume.
+
+{{< img src="../RichardsComponentTransport_Padilla.png" title="Comparison between numerical (lines) and experimental (squares) results for cases 'NaCl1' and 'NaCl6' from Padilla et al. (1999).">}}
+
+It can be seen, that with decreasing saturation, breakthrough curves exhibit stronger dispersion through the decreased angle of the breakthrough curve. Both simulation results follow the experimental observations closely; deviations, especially in the unsaturated case, can be attributed to known tailing effects from secondary porosity.
+
+Here are digitized data from the [NaCl1](../../../../../Tests/Data/Parabolic/RichardsComponentTransport/Padilla_NaCl1_poreVolume-conc.csv) and [NaCl6](../../../../../Tests/Data/Parabolic/RichardsComponentTransport/Padilla_NaCl6_poreVolume-conc.csv) experiments and a [ParaView state file](../../../../../Tests/Data/Parabolic/RichardsComponentTransport/Padilla_state.pvsm) for comparison.
+
+
+## References
+
+Padilla, I. Y., T.-C. J. Yeh, and M. H. Conklin (1999), The effect of water content on solute transport in unsaturated porous media, Water Resour. Res., 35(11), 3303–3313, doi:10.1029/1999WR900171.
+