diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index d214974d3f142273449deaf60706d42d615a0fda..d995718043a40a508a712c0fac17c2836757987c 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -465,7 +465,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                 ProcessLib::ComponentTransport::createComponentTransportProcess(
                     *_mesh_vec[0], std::move(jacobian_assembler),
                     _process_variables, _parameters, integration_order,
-                    process_config);
+                    process_config, _mesh_vec, output_directory);
         }
         else
 #endif
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 0558f4f3941c717d9d74f0e6ff1a1a78832fcbc7..a5f63eacca9734fc9b45bf03cf6f34b3398c3960 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -336,6 +336,70 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
+    Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
+                            double const t,
+                            std::vector<double> const& local_x) const override
+    {
+        // eval shape matrices at given point
+        auto const shape_matrices = [&]() {
+            using FemType =
+                NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+            FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+                &_element));
+
+            typename ShapeMatricesType::ShapeMatrices shape_matrices(
+                ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+            // Note: Axial symmetry is set to false here, because we only need
+            // dNdx here, which is not affected by axial symmetry.
+            fe.template computeShapeFunctions<NumLib::ShapeMatrixType::DNDX>(
+                pnt_local_coords.getCoords(), shape_matrices, GlobalDim, false);
+
+            return shape_matrices;
+        }();
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+        // local_x contains the local concentration and pressure values
+        NumLib::shapeFunctionInterpolate(
+            local_x, shape_matrices.N,
+            vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::C)],
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::p)]);
+
+        auto const K =
+            _process_data.porous_media_properties
+                .getIntrinsicPermeability(t, pos)
+                .getValue(t, pos, 0.0,
+                          vars[static_cast<int>(
+                              MaterialLib::Fluid::PropertyVariableType::C)]);
+
+        auto const mu = _process_data.fluid_properties->getValue(
+            MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+        GlobalDimMatrixType const K_over_mu = K / mu;
+
+        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+            &local_x[local_x.size()/2], ShapeFunction::NPOINTS);
+        GlobalDimVectorType q =
+            -K_over_mu * shape_matrices.dNdx * p_nodal_values;
+
+        if (_process_data.has_gravity)
+        {
+            auto const rho_w =
+                _process_data.fluid_properties->getValue(
+                    MaterialLib::Fluid::FluidPropertyType::Density, vars);
+            auto const b = _process_data.specific_body_force;
+            q += K_over_mu * rho_w * b;
+        }
+        Eigen::Vector3d flux(0.0, 0.0, 0.0);
+        flux.head<GlobalDim>() = q;
+        return flux;
+    }
+
 private:
     MeshLib::Element const& _element;
     ComponentTransportProcessData const& _process_data;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 0d8a1a4eb2d357d8a6ddb0c833f7fd825f9e284a..6ba33a623fad21cce79d4de4a10620de2787d384 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -11,6 +11,8 @@
 
 #include <cassert>
 
+#include "ProcessLib/SurfaceFlux/SurfaceFlux.h"
+#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
 
 namespace ProcessLib
@@ -27,12 +29,14 @@ ComponentTransportProcess::ComponentTransportProcess(
     ComponentTransportProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
     NumLib::NamedFunctionCaller&& named_function_caller,
-    bool const use_monolithic_scheme)
+    bool const use_monolithic_scheme,
+    std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux)
     : Process(mesh, std::move(jacobian_assembler), parameters,
               integration_order, std::move(process_variables),
               std::move(secondary_variables), std::move(named_function_caller),
               use_monolithic_scheme),
-      _process_data(std::move(process_data))
+      _process_data(std::move(process_data)),
+      _surfaceflux(std::move(surfaceflux))
 {
 }
 
@@ -88,6 +92,46 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, _coupled_solutions);
 }
 
+Eigen::Vector3d ComponentTransportProcess::getFlux(std::size_t const element_id,
+                                                   MathLib::Point3d const& p,
+                                                   double const t,
+                                                   GlobalVector const& x) const
+{
+    std::vector<GlobalIndexType> indices_cache;
+    auto const r_c_indices = NumLib::getRowColumnIndices(
+        element_id, *_local_to_global_index_map, indices_cache);
+    std::vector<double> local_x(x.get(r_c_indices.rows));
+    return _local_assemblers[element_id]->getFlux(p, t, local_x);
+}
+
+void ComponentTransportProcess::postTimestepConcreteProcess(
+    GlobalVector const& x,
+    const double t,
+    const double /*delta_t*/,
+    int const process_id)
+{
+    // For the monolithic scheme, process_id is always zero.
+    if (_use_monolithic_scheme && process_id != 0)
+    {
+        OGS_FATAL(
+            "The condition of process_id = 0 must be satisfied for "
+            "monolithic ComponentTransportProcess, which is a single process.");
+    }
+    if (!_use_monolithic_scheme && process_id != 1)
+    {
+        DBUG(
+            "This is the transport part of the staggered "
+            "ComponentTransportProcess.");
+        return;
+    }
+    if (!_surfaceflux)  // computing the surfaceflux is optional
+    {
+        return;
+    }
+    _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh);
+    _surfaceflux->save(t);
+}
+
 }  // namespace ComponentTransport
 }  // namespace ProcessLib
 
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 03dbc62af4e5fd04c97735febb1cc9c5e9dc4085..78e6905f71974ff10454e874cfd2bf671a9a9439 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -16,6 +16,8 @@
 
 namespace ProcessLib
 {
+struct SurfaceFluxData;
+
 namespace ComponentTransport
 {
 /**
@@ -97,7 +99,8 @@ public:
         ComponentTransportProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,
         NumLib::NamedFunctionCaller&& named_function_caller,
-        bool const use_monolithic_scheme);
+        bool const use_monolithic_scheme,
+        std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux);
 
     //! \name ODESystem interface
     //! @{
@@ -105,6 +108,15 @@ public:
     bool isLinear() const override { return false; }
     //! @}
 
+    Eigen::Vector3d getFlux(std::size_t const element_id,
+                            MathLib::Point3d const& p, double const t,
+                            GlobalVector const& x) const override;
+
+    void postTimestepConcreteProcess(GlobalVector const& x,
+                                     const double t,
+                                     const double delta_t,
+                                     int const process_id) override;
+
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -124,6 +136,8 @@ private:
 
     std::vector<std::unique_ptr<ComponentTransportLocalAssemblerInterface>>
         _local_assemblers;
+
+    std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux;
 };
 
 }  // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index c7c73b9ff219087d1cd00a932d0332e83f052dd7..a5c3f2663c0f7e06b704b8457fad750925a297d5 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -12,6 +12,9 @@
 #include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
 #include "MaterialLib/PorousMedium/CreatePorousMediaProperties.h"
 
+#include "MeshLib/IO/readMeshFromFile.h"
+
+#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Parameter/ConstantParameter.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
@@ -29,7 +32,9 @@ std::unique_ptr<Process> createComponentTransportProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
+    std::string const& output_directory)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "ComponentTransport");
@@ -165,11 +170,21 @@ std::unique_ptr<Process> createComponentTransportProcess(
     ProcessLib::createSecondaryVariables(config, secondary_variables,
                                          named_function_caller);
 
+    std::unique_ptr<ProcessLib::SurfaceFluxData> surfaceflux;
+    auto surfaceflux_config =
+        //! \ogs_file_param{prj__processes__process__calculatesurfaceflux}
+        config.getConfigSubtreeOptional("calculatesurfaceflux");
+    if (surfaceflux_config)
+    {
+        surfaceflux = ProcessLib::SurfaceFluxData::createSurfaceFluxData(
+            *surfaceflux_config, meshes, output_directory);
+    }
+
     return std::make_unique<ComponentTransportProcess>(
         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),
-        use_monolithic_scheme);
+        use_monolithic_scheme, std::move(surfaceflux));
 }
 
 }  // namespace ComponentTransport
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
index 98071011f2634d24033efe25994218ae9489f419..8249620bb7ab71c4e57902d67f9428c54bdf0a74 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
@@ -22,7 +22,9 @@ std::unique_ptr<Process> createComponentTransportProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
+    std::string const& output_directory);
 
 }  // namespace ComponentTransport
 }  // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index c0075687e6764c9ae877d6894649ad245d4c6f4e..564882f4a5231767a11c56982ceaf37a105dee8e 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -193,6 +193,18 @@ AddTest(
     VIS DiffusionAndStorageAndAdvectionAndDispersionHalf_pcs_0_ts_672_t_900.000000.vtu
 )
 
+AddTest(
+    NAME 3D_ComponentTransport_surfaceflux
+    PATH Parabolic/ComponentTransport/SimpleSynthetics
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS surfaceflux_component-transport_cube_1e3.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    cube_1x1x1_hex_1e3_complete_surface.vtu flux_1e3_t_1.000000.vtu flux flux 1e-10 1e-16
+)
+
 AddTest(
     NAME LARGE_2D_ComponentTransport_Goswami
     PATH Parabolic/ComponentTransport/goswami
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index dc179be6c366fcf733e17870858f4cc5aaddd98b..a80ff47d368d38e7ac37cfc6537e2d0ea7c1a6a8 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -13,11 +13,10 @@
 
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
-
+#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
 
 #include "HTMaterialProperties.h"
-
 #include "MonolithicHTFEM.h"
 #include "StaggeredHTFEM.h"
 
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index b9f74a6b9bd4c892c25dee614f54fa8a453216cd..b4c9038748d0d619157723e9864c14eb0c0bf49e 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -11,8 +11,6 @@
 
 #include <array>
 
-#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
-#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
 #include "ProcessLib/Process.h"
 
 namespace NumLib
@@ -22,6 +20,8 @@ class LocalToGlobalIndexMap;
 
 namespace ProcessLib
 {
+struct SurfaceFluxData;
+
 namespace HT
 {
 class HTLocalAssemblerInterface;
diff --git a/ProcessLib/SurfaceFlux/SurfaceFlux.h b/ProcessLib/SurfaceFlux/SurfaceFlux.h
index 4101b6ba6ff29673ffac9be6358424f7f733952e..1c9eaf74887ae35684ad9058f74be54e329a85e4 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFlux.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFlux.h
@@ -23,8 +23,8 @@ public:
     /// variable has.
     /// @param integration_order Integration order used in local assembly.
     SurfaceFlux(MeshLib::Mesh& boundary_mesh,
-                         std::size_t bulk_property_number_of_components,
-                         unsigned const integration_order);
+                std::size_t bulk_property_number_of_components,
+                unsigned const integration_order);
 
     /// Executes for each element of the mesh the local intergration procedure.
     /// @param x The global solution the intergration values will be fetched of.
diff --git a/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3.vtu b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..5f491d77c01fe3246feba3443d5548791e934e82
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:61a70995845833f16ba3b2edf0312eca29868be20f2cc7015df06d1f35eed669
+size 22546
diff --git a/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_complete_surface.vtu b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_complete_surface.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8a2407e5f3183a2091ee355e342e9c8ab505e457
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_complete_surface.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bcb45ed892b69e6ddf23038d5d6dbe292a2f95448b97e6401336bbdb4f8982e5
+size 79181
diff --git a/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_left.vtu b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_left.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..b9be1a60800686b8a6c073613bd5169f6f69ccad
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_left.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:074db4f230b18b0fc6300e7b49669b55cd2e06b948b1344302305a59261de79d
+size 12897
diff --git a/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_right.vtu b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_right.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8c9a97647805e21f6efe641cec909daae8b5f394
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/cube_1x1x1_hex_1e3_right.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8427200cf2b82af3762994c1e259839aa5b502a8d5b8a7dce60919f70fa04e41
+size 12898
diff --git a/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/surfaceflux_component-transport_cube_1e3.prj b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/surfaceflux_component-transport_cube_1e3.prj
new file mode 100644
index 0000000000000000000000000000000000000000..e0a53624697ec24f631eacd9d1d32af6b456ab7d
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/SimpleSynthetics/surfaceflux_component-transport_cube_1e3.prj
@@ -0,0 +1,240 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>cube_1x1x1_hex_1e3.vtu</mesh>
+        <mesh>cube_1x1x1_hex_1e3_left.vtu</mesh>
+        <mesh>cube_1x1x1_hex_1e3_right.vtu</mesh>
+        <mesh>cube_1x1x1_hex_1e3_complete_surface.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>hc</name>
+            <type>ComponentTransport</type>
+            <integration_order>2</integration_order>
+            <process_variables>
+                <concentration>concentration</concentration>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <fluid>
+                <density>
+                    <type>ConcentrationDependent</type>
+                    <reference_density>1</reference_density>
+                    <reference_concentration>0</reference_concentration>
+                    <fluid_density_difference_ratio>0.0</fluid_density_difference_ratio>
+                </density>
+                <viscosity>
+                    <type>Constant</type>
+                    <value>1.0e-3</value>
+                </viscosity>
+            </fluid>
+            <porous_medium>
+                <porous_medium id="0">
+                    <permeability>
+                        <type>Constant</type>
+                        <permeability_tensor_entries>kappa1</permeability_tensor_entries>
+                    </permeability>
+                    <porosity>
+                        <type>Constant</type>
+                        <porosity_parameter>constant_porosity_parameter</porosity_parameter>
+                    </porosity>
+                    <storage>
+                        <type>Constant</type>
+                        <value>0</value>
+                    </storage>
+                </porous_medium>
+            </porous_medium>
+            <fluid_reference_density>rho_fluid</fluid_reference_density>
+            <solute_dispersivity_longitudinal>beta_l</solute_dispersivity_longitudinal>
+            <solute_dispersivity_transverse>beta_t</solute_dispersivity_transverse>
+            <molecular_diffusion_coefficient>Dm</molecular_diffusion_coefficient>
+            <retardation_factor>retardation</retardation_factor>
+            <decay_rate>decay</decay_rate>
+            <specific_body_force>0 0 0</specific_body_force>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="darcy_velocity" output_name="darcy_velocity"/>
+            </secondary_variables>
+            <calculatesurfaceflux>
+                <mesh>cube_1x1x1_hex_1e3_complete_surface</mesh>
+                <property_name>flux</property_name>
+                <output_mesh>flux_1e3.vtu</output_mesh>
+            </calculatesurfaceflux>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="hc">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1e-1</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <output>
+                    <variables>
+                        <variable>concentration</variable>
+                        <variable>pressure</variable>
+                        <variable>darcy_velocity</variable>
+                    </variables>
+                </output>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0.0</t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>DiffusionOnly_surfaceflux</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>kappa1</name>
+            <type>Constant</type>
+            <values>1.239e-7 0 0 0 1.239e-7 0 0 0 1.239e-7</values>
+        </parameter>
+        <parameter>
+            <name>constant_porosity_parameter</name>
+            <type>Constant</type>
+            <value>0.2</value>
+        </parameter>
+        <parameter>
+            <name>rho_fluid</name>
+            <type>Constant</type>
+            <value>1000</value>
+        </parameter>
+        <parameter>
+            <name>Dm</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>retardation</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>decay</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>beta_l</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>beta_t</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>c0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p_Dirichlet_left</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>p_Dirichlet_right</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>c_Dirichlet_left</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>c_Dirichlet_right</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>concentration</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>c0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>cube_1x1x1_hex_1e3_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_Dirichlet_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>cube_1x1x1_hex_1e3_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>c_Dirichlet_right</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>cube_1x1x1_hex_1e3_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>cube_1x1x1_hex_1e3_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet_right</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>2</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 20000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>ILUT</precon_type>
+                <max_iteration_step>20000</max_iteration_step>
+                <error_tolerance>1e-8</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>hc</prefix>
+                <parameters>-hc_ksp_type bcgs -hc_pc_type bjacobi -hc_ksp_rtol 1e-8 -hc_ksp_max_it 20000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>