diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index c3b5e62c9f39832eb9399a1d78ed9e4639a87f59..e449fe278eefcb456dfeff6581b8fc14186a4201 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -538,6 +538,38 @@ public:
         }
     }
 
+    GlobalDimMatrixType getHydrodynamicDispersion(
+        GlobalDimMatrixType const& I,
+        GlobalDimMatrixType const& pore_diffusion_coefficient,
+        GlobalDimVectorType const& velocity,
+        double const porosity,
+        double const solute_dispersivity_transverse,
+        double const solute_dispersivity_longitudinal)
+    {
+        double const velocity_magnitude = velocity.norm();
+
+        if (velocity_magnitude == 0.0)
+        {
+            return porosity * pore_diffusion_coefficient;
+        }
+
+        GlobalDimMatrixType const D_c = GlobalDimMatrixType(
+            porosity * pore_diffusion_coefficient +
+            solute_dispersivity_transverse * velocity_magnitude * I +
+            (solute_dispersivity_longitudinal -
+             solute_dispersivity_transverse) /
+                velocity_magnitude * velocity * velocity.transpose());
+
+        if (!_process_data.stabilizer)
+        {
+            return D_c;
+        }
+
+        return D_c + _process_data.stabilizer->getExtraDiffusionCoefficient(
+                         _element.getID(), 1.0, velocity_magnitude) *
+                         I;
+    }
+
     void assembleBlockMatrices(
         int const component_id, double const t, double const dt,
         Eigen::Ref<const NodalVectorType> const& C_nodal_values,
@@ -655,24 +687,16 @@ public:
                         vars, MaterialPropertyLib::Variable::concentration, pos,
                         t, dt);
 
-            double const velocity_magnitude = velocity.norm();
             GlobalDimMatrixType const hydrodynamic_dispersion =
-                velocity_magnitude != 0.0
-                    ? GlobalDimMatrixType(porosity *
-                                              pore_diffusion_coefficient +
-                                          solute_dispersivity_transverse *
-                                              velocity_magnitude * I +
-                                          (solute_dispersivity_longitudinal -
-                                           solute_dispersivity_transverse) /
-                                              velocity_magnitude * velocity *
-                                              velocity.transpose())
-                    : GlobalDimMatrixType(porosity *
-                                              pore_diffusion_coefficient +
-                                          solute_dispersivity_transverse *
-                                              velocity_magnitude * I);
+                getHydrodynamicDispersion(I, pore_diffusion_coefficient,
+                                          velocity, porosity,
+                                          solute_dispersivity_transverse,
+                                          solute_dispersivity_longitudinal);
+
             const double R_times_phi(retardation_factor * porosity);
             GlobalDimVectorType const mass_density_flow = velocity * density;
             auto const N_t_N = (N.transpose() * N).eval();
+
             if (_process_data.non_advective_form)
             {
                 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
@@ -1020,21 +1044,11 @@ public:
                                           (dNdx * local_p - density * b))
                     : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
 
-            double const velocity_magnitude = velocity.norm();
             GlobalDimMatrixType const hydrodynamic_dispersion =
-                velocity_magnitude != 0.0
-                    ? GlobalDimMatrixType(porosity *
-                                              pore_diffusion_coefficient +
-                                          solute_dispersivity_transverse *
-                                              velocity_magnitude * I +
-                                          (solute_dispersivity_longitudinal -
-                                           solute_dispersivity_transverse) /
-                                              velocity_magnitude * velocity *
-                                              velocity.transpose())
-                    : GlobalDimMatrixType(porosity *
-                                              pore_diffusion_coefficient +
-                                          solute_dispersivity_transverse *
-                                              velocity_magnitude * I);
+                getHydrodynamicDispersion(I, pore_diffusion_coefficient,
+                                          velocity, porosity,
+                                          solute_dispersivity_transverse,
+                                          solute_dispersivity_longitudinal);
 
             double const R_times_phi = retardation_factor * porosity;
             auto const N_t_N = (N.transpose() * N).eval();
@@ -1053,7 +1067,6 @@ public:
             local_M.noalias() += N_t_N * (R_times_phi * density * w);
 
             // coupling term
-
             if (_process_data.non_advective_form)
             {
                 double dot_p_int_pt = 0.0;
@@ -1317,14 +1330,8 @@ public:
                     ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
                     : GlobalDimVectorType(-k / mu * dNdx * p);
 
-            double const q_magnitude = q.norm();
-            // hydrodymanic dispersion
             GlobalDimMatrixType const D =
-                q_magnitude != 0.0
-                    ? GlobalDimMatrixType(phi * Dp + alpha_T * q_magnitude * I +
-                                          (alpha_L - alpha_T) / q_magnitude *
-                                              q * q.transpose())
-                    : GlobalDimMatrixType(phi * Dp);
+                getHydrodynamicDispersion(I, Dp, q, phi, alpha_T, alpha_L);
 
             // matrix assembly
             local_Jac.noalias() +=
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
index b022a186420b77e24a412bc3957eb3b6869faecd..33d3f91cd660c098f71146bf3de3d4b51b136e1e 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
@@ -16,6 +16,7 @@
 #include "LookupTable.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/NumericalStability/NumericalStabilization.h"
 #include "ParameterLib/Parameter.h"
 
 namespace MaterialPropertyLib
@@ -69,6 +70,8 @@ struct ComponentTransportProcessData
     ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface;
     std::unique_ptr<LookupTable> lookup_table;
 
+    std::unique_ptr<NumLib::NumericalStabilization> stabilizer;
+
     static const int hydraulic_process_id = 0;
     // TODO (renchao-lu): This variable is used in the calculation of the
     // fluid's density and flux, indicating the transport process id. For now it
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index 873081551d4d9730621ca3e0ece30bcb637b6d47..e525ac72e62e10c5226ec9046e05e52f055aa1f0 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -236,6 +236,7 @@ std::unique_ptr<Process> createComponentTransportProcess(
     checkMPLProperties(mesh, *media_map);
     DBUG("Media properties verified.");
 
+    auto stabilizer = NumLib::createNumericalStabilization(mesh, config);
     ComponentTransportProcessData process_data{
         std::move(media_map),
         specific_body_force,
@@ -244,7 +245,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
         temperature,
         chemically_induced_porosity_change,
         chemical_solver_interface.get(),
-        std::move(lookup_table)};
+        std::move(lookup_table),
+        std::move(stabilizer)};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index 20ed145c85cd8f662440588f9ea1437019dd00b5..af60513ad11bbe5dd5963e91f2fc9437b38d8db5 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -845,3 +845,21 @@ if (OGS_USE_MPI)
     OgsTest(WRAPPER mpirun -np 1 PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/EquilibriumPhase/calcitePorosityChange.prj RUNTIME 25)
     OgsTest(WRAPPER mpirun -np 2 PROJECTFILE Parabolic/ComponentTransport/ReactiveTransport/SurfaceComplexation/ParallelTest/RadionuclideSorption.prj RUNTIME 60)
 endif()
+
+AddTest(
+    NAME ComponentTransport_ClassicalTransportExample
+    PATH Parabolic/ComponentTransport/ClassicalTransportExample
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS classical_transport_example.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    RUNTIME 1
+    DIFF_DATA
+    classical_transport_example_t_4800.00.vtu classical_transport_example_t_4800.00.vtu C C 1.e-9 1.0e-12
+    classical_transport_example_t_4800.00.vtu classical_transport_example_t_4800.00.vtu pressure pressure 1.e-9 1.0e-12
+    classical_transport_example_t_4800.00.vtu classical_transport_example_t_4800.00.vtu velocity velocity 1.e-12 1.0e-12
+    classical_transport_example_t_7200.00.vtu classical_transport_example_t_7200.00.vtu C C 1.e-9 1.0e-12
+    classical_transport_example_t_7200.00.vtu classical_transport_example_t_7200.00.vtu pressure pressure 1.e-9 1.0e-12
+    classical_transport_example_t_7200.00.vtu classical_transport_example_t_7200.00.vtu velocity velocity 1.e-12 1.0e-12
+)
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example.prj b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example.prj
new file mode 100644
index 0000000000000000000000000000000000000000..e53048e695c8c6abaa8afe163d6d154d07b87358
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example.prj
@@ -0,0 +1,221 @@
+<?xml version='1.0' encoding='ISO-8859-1'?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>mesh_2D.vtu</mesh>
+        <mesh>geometry_right.vtu</mesh>
+        <mesh>geometry_left.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>HC</name>
+            <type>ComponentTransport</type>
+            <integration_order>4</integration_order>
+            <process_variables>
+                <concentration>C</concentration>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="darcy_velocity" output_name="darcy_velocity"/>
+            </secondary_variables>
+            <specific_body_force>0 0</specific_body_force>
+            <numerical_stabilization>
+                <type>IsotropicDiffusion</type>
+                <tuning_parameter>0.15</tuning_parameter>
+                <cutoff_velocity>0.e-0</cutoff_velocity>
+            </numerical_stabilization>
+        </process>
+    </processes>
+    <media>
+        <medium id="0">
+            <phases>
+                <phase>
+                    <type>AqueousLiquid</type>
+                    <components>
+                        <component>
+                            <name>C</name>
+                            <properties>
+                                <property>
+                                    <name>pore_diffusion</name>
+                                    <type>Constant</type>
+                                    <value>1.e-9</value>
+                                </property>
+                                <property>
+                                    <name>retardation_factor</name>
+                                    <type>Constant</type>
+                                    <value>1.0</value>
+                                </property>
+                                <property>
+                                    <name>decay_rate</name>
+                                    <type>Parameter</type>
+                                    <parameter_name>decay</parameter_name>
+                                </property>
+                            </properties>
+                        </component>
+                    </components>
+                    <properties>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value>1</value>
+                        </property>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value>1.0</value>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>permeability</name>
+                    <type>Constant</type>
+                    <value>1.e-9</value>
+                </property>
+                <property>
+                    <name>porosity</name>
+                    <type>Constant</type>
+                    <value>1.0</value>
+                </property>
+                <property>
+                    <name>longitudinal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0</value>
+                </property>
+                <property>
+                    <name>transversal_dispersivity</name>
+                    <type>Constant</type>
+                    <value>0.0</value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <parameters>
+        <parameter>
+            <name>decay</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>C0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>C1</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p1</name>
+            <type>Constant</type>
+            <value>8.0e+4</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>geometry_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p1</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>geometry_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>C</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>C0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>geometry_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>C0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>geometry_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>C1</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <time_loop>
+        <processes>
+            <process ref="HC">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltols>1e-12 1e-10</reltols>
+                    <abstols>1e-9 1e-6</abstols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial> 0.0 </t_initial>
+                    <t_end> 7200 </t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>200</repeat>
+                            <delta_t>18</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>classical_transport_example</prefix>
+            <suffix>_t_{:0.2ftime}</suffix>
+            <timesteps>
+                <pair>
+                    <repeat>1000</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>C</variable>
+                <variable>pressure</variable>
+                <variable>velocity</variable>
+            </variables>
+            <fixed_output_times>
+                3600 4800
+            </fixed_output_times>
+        </output>
+    </time_loop>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>50</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <scaling>true</scaling>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_4800.00.vtu b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_4800.00.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..baf082daf0d97b593d9f51133de705963465cb7d
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_4800.00.vtu
@@ -0,0 +1,32 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="22" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+    </FieldData>
+    <Piece NumberOfPoints="30"                   NumberOfCells="14"                  >
+      <PointData>
+        <DataArray type="Float64" Name="C" format="appended" RangeMin="-0.0022197086164"     RangeMax="1.0020435168"         offset="84"                  />
+        <DataArray type="Float64" Name="MolarFlowRate" format="appended" RangeMin="0"                    RangeMax="0"                    offset="380"                 />
+        <DataArray type="Float64" Name="darcy_velocity" NumberOfComponents="2" format="appended" RangeMin="0.0001"               RangeMax="0.0001"               offset="440"                 />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="0"                    RangeMax="80000"                offset="968"                 />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1244"                />
+        <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1304"                />
+        <DataArray type="Float64" Name="velocity" NumberOfComponents="2" format="appended" RangeMin="0.0001"               RangeMax="0.0001"               offset="1364"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.80006249756"        offset="1540"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1888"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="2056"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="2160"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABYAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9I1NDS11E03tzCwMDA3AmIANB4EmQ==AQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAAvAAAAAAAAAA=eF5jYACBD/YMWMEH+3vMm9Y4/Htvb9G8TGA70wf72e8jAy79f2+fW2Y8Ne35e/v8WTEn4jg+2DefdonQXPzevqfrgUsj+xv7VVe33Z+i88R+/fS6G1O2XLVPUF38Q691v72t9kKGpe7T7T/a2f0V0E/aH+NzTNNaKHH/2Q1HwfSLHfZg8cwlc8Hqmtvmg/UZz6sBm3OuYyvY3N3T74PtOT4NYu/yuRB37IW66+JbiDtzoO5+DfUHAKDmbf8=AQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAADAAAAAAAAAA=eF5jYBhZAAAA8AABAQAAAAAAAAAAgAAAAAAAAOABAAAAAAAAaQEAAAAAAAA=eF6r9ZF5/chMyj7qS3ZtwWJHmzw7CH9GgXmI1crfu4uhfO9Uu4nJFn93p0DVbwrZdPvFVR+bcCh/f4v9TVddzT3vvaDqNx/N7rIw35PrAeEnW6far0wR39PkBuGft/7q8bZEf89iZwh/gj57bMqyJ9ZBjhD+Emu9F+4NKnu+OED4pqdblfXma9n02UP4OdOfdPgtNrF5CHXf7fsHgtwmCuxpgspfvKzQWb9CwiYCyq/JUPXRm6i0hwHK3zCfk+mQnaTNJqj+GheOytR313dvg/HdGgO89f7v5oGq73rO/9N9lbBNApRvcHjdGtFUpT0dUP6fDf/N1bSYbB5D9Ufd69q+9duH3Z1Q+RkhWcfmdWrZ/IH6p+JoYmB0gKJNANS/MxKEbKte/rdeAA0Pr9Jf2Wr6UjYt0PCykZjXvHqn9p4SaHjOYZnyJOOk5Z6X0PDOk3ydYldtvscKGh87prht/ZDNtQcAZD/S8g==AQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAArgAAAAAAAAA=eF5jYACCjs8ODGjgPwi0f3a4dndb7m2lTw6tu3Nvb9v7wWG6mecklY3vHC5/BYo/e+1gulXFc5LMSwfLskkqnkHPHJpNQeY9dnj+BKT+gcNU9e1AdbcdyttA8tcc/J+CxC84dHWA+MccNnSD6G0O+owQOo8RIq7WCFG3mRGiz6RpG9gc+0aIuRaHQO577GAtDLGXMRXijunbIO7y54C4U2MtxN0x1yD+AACMOXGJAQAAAAAAAAAAgAAAAAAAADgAAAAAAAAADAAAAAAAAAA=eF5jYCAPAAAAOAABAQAAAAAAAAAAgAAAAAAAAHAAAAAAAAAADAAAAAAAAAA=eF5jYKAvAAAAcAABAQAAAAAAAAAAgAAAAAAAAOAAAAAAAAAAYgAAAAAAAAA=eF7L95F5/chMyp4BBBr0bORR+Ep7tnki8QXM92i5I/EZdPZcdUXiO4juyXBElv+w2wyFL2qz2QHVvjg7ZPkTu5Pskfn8NotQ+Ox7BJD5Caw2P5D1JzDsyUcxb481ACgROhM=AQAAAAAAAAAAgAAAAAAAANACAAAAAAAA4wAAAAAAAAA=eF5jYMAOZs0EgZf2uMSrRda5P6xqwZBHF3ec7ih7xXEthjrOmSDxvRjiHzU+isZ/PIoh/rQNpP4shnj2J49JKp6XMMTvyoPMuYrpfn+Q+29iiNtPAZl/F0P8w5uSAPbSBxjiQfyeQHsfYYhvMgQaI/sEQ1w4BOSepxjiE8qOhtccfY4h/qIRIo4enmZZEHPQxVuCIPaii1+2grgTXVxRuRTsL3Txu2JO4HBAF9/2DBJu6OLzFkLCGV28O9oTHC/o4mt0nMDxiC7uuggS7+ji6/WdwOkEXVzAxAmcrtDFASQHob4=AQAAAAAAAAAAgAAAAAAAAMABAAAAAAAAWwAAAAAAAAA=eF5dyrUNwAAAA8EwM+P+a6bIu7Gbk6wPgn8Jvhjbn+Jjnf4Mb+v053hZp7/A0zr9JR7W6a9wt05/jZt1+htcrdPf4mKd/g5n6/T3OFmnf8DROv0hRtZ9+KYDYQ==AQAAAAAAAAAAgAAAAAAAAHAAAAAAAAAALQAAAAAAAAA=eF5jYYAADijNA6UFoLQIlJaA0jJQWgFKq0BpDSitA6UNoLQJlLaA0gBGcAGlAQAAAAAAAAAAgAAAAAAAAA4AAAAAAAAACwAAAAAAAAA=eF7j5EQGAAO/AH8=
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_7200.00.vtu b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_7200.00.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..583de7471307a6f323e81c8448a1b15c9cad9f3e
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example_t_7200.00.vtu
@@ -0,0 +1,32 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="22" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
+    </FieldData>
+    <Piece NumberOfPoints="30"                   NumberOfCells="14"                  >
+      <PointData>
+        <DataArray type="Float64" Name="C" format="appended" RangeMin="0"                    RangeMax="1.0060206808"         offset="84"                  />
+        <DataArray type="Float64" Name="MolarFlowRate" format="appended" RangeMin="0"                    RangeMax="0"                    offset="364"                 />
+        <DataArray type="Float64" Name="darcy_velocity" NumberOfComponents="2" format="appended" RangeMin="0.0001"               RangeMax="0.0001"               offset="424"                 />
+        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="0"                    RangeMax="80000"                offset="952"                 />
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1228"                />
+        <DataArray type="Float64" Name="porosity_avg" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1288"                />
+        <DataArray type="Float64" Name="velocity" NumberOfComponents="2" format="appended" RangeMin="0.0001"               RangeMax="0.0001"               offset="1348"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.80006249756"        offset="1524"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1872"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="2040"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="2144"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAABYAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9I1NDS11E03tzCwMDA3AmIANB4EmQ==AQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAAsQAAAAAAAAA=eF77/x8E3tszoIG/UPGHm1MY3wLprKVrDoPovO+rrfoYPthzCHlNZ/r73t7w12xtRqYP9vX8xvNs/7y3/3js1ZxpLB/sl857/TL74nt73/N6KislPtgvefLlqm3RO/up3u1n7k54Y+93o0e0VfmxfXDdqnM2Ox/Y9/tD6FWdEPG5t9rA6hKOQvQFHYGYE7MSYu6RQxB7GoUg9q75CXHHFqi7bKDunAh19zWoPwDvLIBYAQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAADAAAAAAAAAA=eF5jYBhZAAAA8AABAQAAAAAAAAAAgAAAAAAAAOABAAAAAAAAaQEAAAAAAAA=eF6r9ZF5/chMyj7qS3ZtwWJHmzw7CH9GgXmI1crfu4uhfO9Uu4nJFn93p0DVbwrZdPvFVR+bcCh/f4v9TVddzT3vvaDqNx/N7rIw35PrAeEnW6far0wR39PkBuGft/7q8bZEf89iZwh/gj57bMqyJ9ZBjhD+Emu9F+4NKnu+OED4pqdblfXma9n02UP4OdOfdPgtNrF5CHXf7fsHgtwmCuxpgspfvKzQWb9CwiYCyq/JUPXRm6i0hwHK3zCfk+mQnaTNJqj+GheOytR313dvg/HdGgO89f7v5oGq73rO/9N9lbBNApRvcHjdGtFUpT0dUP6fDf/N1bSYbB5D9Ufd69q+9duH3Z1Q+RkhWcfmdWrZ/IH6p+JoYmB0gKJNANS/MxKEbKte/rdeAA0Pr9Jf2Wr6UjYt0PCykZjXvHqn9p4SaHjOYZnyJOOk5Z6X0PDOk3ydYldtvscKGh87prht/ZDNtQcAZD/S8g==AQAAAAAAAAAAgAAAAAAAAPAAAAAAAAAArgAAAAAAAAA=eF5jYACCjs8ODGjgPwi0f3a4dndb7m2lTw6tu3Nvb9v7wWG6mecklY3vHC5/BYo/e+1gulXFc5LMSwfLskkqnkHPHJpNQeY9dnj+BKT+gcNU9e1AdbcdyttA8tcc/J+CxC84dHWA+MccNnSD6G0O+owQOo8RIq7WCFG3mRGiz6RpG9gc+0aIuRaHQO577GAtDLGXMRXijunbIO7y54C4U2MtxN0x1yD+AACMOXGJAQAAAAAAAAAAgAAAAAAAADgAAAAAAAAADAAAAAAAAAA=eF5jYCAPAAAAOAABAQAAAAAAAAAAgAAAAAAAAHAAAAAAAAAADAAAAAAAAAA=eF5jYKAvAAAAcAABAQAAAAAAAAAAgAAAAAAAAOAAAAAAAAAAYgAAAAAAAAA=eF7L95F5/chMyp4BBBr0bORR+Ep7tnki8QXM92i5I/EZdPZcdUXiO4juyXBElv+w2wyFL2qz2QHVvjg7ZPkTu5Pskfn8NotQ+Ox7BJD5Caw2P5D1JzDsyUcxb481ACgROhM=AQAAAAAAAAAAgAAAAAAAANACAAAAAAAA4wAAAAAAAAA=eF5jYMAOZs0EgZf2uMSrRda5P6xqwZBHF3ec7ih7xXEthjrOmSDxvRjiHzU+isZ/PIoh/rQNpP4shnj2J49JKp6XMMTvyoPMuYrpfn+Q+29iiNtPAZl/F0P8w5uSAPbSBxjiQfyeQHsfYYhvMgQaI/sEQ1w4BOSepxjiE8qOhtccfY4h/qIRIo4enmZZEHPQxVuCIPaii1+2grgTXVxRuRTsL3Txu2JO4HBAF9/2DBJu6OLzFkLCGV28O9oTHC/o4mt0nMDxiC7uuggS7+ji6/WdwOkEXVzAxAmcrtDFASQHob4=AQAAAAAAAAAAgAAAAAAAAMABAAAAAAAAWwAAAAAAAAA=eF5dyrUNwAAAA8EwM+P+a6bIu7Gbk6wPgn8Jvhjbn+Jjnf4Mb+v053hZp7/A0zr9JR7W6a9wt05/jZt1+htcrdPf4mKd/g5n6/T3OFmnf8DROv0hRtZ9+KYDYQ==AQAAAAAAAAAAgAAAAAAAAHAAAAAAAAAALQAAAAAAAAA=eF5jYYAADijNA6UFoLQIlJaA0jJQWgFKq0BpDSitA6UNoLQJlLaA0gBGcAGlAQAAAAAAAAAAgAAAAAAAAA4AAAAAAAAACwAAAAAAAAA=eF7j5EQGAAO/AH8=
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_left.vtu b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_left.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..1ed1022cb55c7192e4fa9a0e929ff23cfaff1af7
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_left.vtu
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="2"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0"                    RangeMax="3"                    offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="1.8446744074e+19"     RangeMax="1.8446744074e+19"     offset="32"                  />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.01"                 offset="56"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="132"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="164"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="188"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _EAAAAAAAAAADAAAAAAAAAAAAAAAAAAAACAAAAAAAAAD//////////w==MAAAAAAAAAAAAAAAAAAAAHsUrkfheoQ/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA=EAAAAAAAAAABAAAAAAAAAAAAAAAAAAAACAAAAAAAAAACAAAAAAAAAA==AQAAAAAAAAAD
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_right.vtu b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_right.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..509dd9aeee392e83f32afe22073b2587fc59dba5
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/geometry_right.vtu
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="2"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="1"                    RangeMax="2"                    offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="13"                   RangeMax="13"                   offset="32"                  />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0.8"                  RangeMax="0.80006249756"        offset="56"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="132"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="164"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="188"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _EAAAAAAAAAACAAAAAAAAAAEAAAAAAAAACAAAAAAAAAANAAAAAAAAAA==MAAAAAAAAACamZmZmZnpP3sUrkfheoQ/AAAAAAAAAACamZmZmZnpPwAAAAAAAAAAAAAAAAAAAAA=EAAAAAAAAAABAAAAAAAAAAAAAAAAAAAACAAAAAAAAAACAAAAAAAAAA==AQAAAAAAAAAD
+  </AppendedData>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/mesh_2D.vtu b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/mesh_2D.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f5236c188379420c0cbff780765472031ad75cef
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ClassicalTransportExample/mesh_2D.vtu
@@ -0,0 +1,23 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="30"                   NumberOfCells="14"                  >
+      <PointData>
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="0"                   />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="0.80006249756"        offset="88"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1060"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="1668"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="1828"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _OAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA==0AIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACamZmZmZnpPwAAAAAAAAAAAAAAAAAAAACamZmZmZnpP3sUrkfheoQ/AAAAAAAAAAAAAAAAAAAAAHsUrkfheoQ/AAAAAAAAAABBl0Ed1EGtPwAAAAAAAAAAAAAAAAAAAAAJmUEd1EG9PwAAAAAAAAAAAAAAAAAAAADxKPEVX/HFPwAAAAAAAAAAAAAAAAAAAADlhkEd1EHNPwAAAAAAAAAAAAAAAAAAAABr8kiSJEnSPwAAAAAAAAAAAAAAAAAAAADdH/EVX/HVPwAAAAAAAAAAAAAAAAAAAACaT5mZmZnZPwAAAAAAAAAAAAAAAAAAAAA/lEEd1EHdPwAAAAAAAAAAAAAAAAAAAADw7HRQB3XgPwAAAAAAAAAAAAAAAAAAAABSD0mSJEniPwAAAAAAAAAAAAAAAAAAAACyMR3UQR3kPwAAAAAAAAAAAAAAAAAAAAATVPEVX/HlPwAAAAAAAAAAAAAAAAAAAACQdsVXfMXnPwAAAAAAAAAAAAAAAAAAAADogcVXfMXnP3sUrkfheoQ/AAAAAAAAAAA2avEVX/HlP3sUrkfheoQ/AAAAAAAAAACEUh3UQR3kP3sUrkfheoQ/AAAAAAAAAADTOkmSJEniP3sUrkfheoQ/AAAAAAAAAAAhI3VQB3XgP3sUrkfheoQ/AAAAAAAAAADdFkId1EHdP3sUrkfheoQ/AAAAAAAAAAC25pmZmZnZP3sUrkfheoQ/AAAAAAAAAACeofEVX/HVP3sUrkfheoQ/AAAAAAAAAACLW0mSJEnSP3sUrkfheoQ/AAAAAAAAAACsLEId1EHNP3sUrkfheoQ/AAAAAAAAAABFovEVX/HFP3sUrkfheoQ/AAAAAAAAAACvL0Id1EG9P3sUrkfheoQ/AAAAAAAAAAAQNEId1EGtP3sUrkfheoQ/AAAAAAAAAAA=wAEAAAAAAAAAAAAAAAAAAAQAAAAAAAAAHQAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAcAAAAAAAAAB0AAAAAAAAABQAAAAAAAAAGAAAAAAAAABsAAAAAAAAAHAAAAAAAAAAGAAAAAAAAAAcAAAAAAAAAGgAAAAAAAAAbAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAZAAAAAAAAABoAAAAAAAAACAAAAAAAAAAJAAAAAAAAABgAAAAAAAAAGQAAAAAAAAAJAAAAAAAAAAoAAAAAAAAAFwAAAAAAAAAYAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAWAAAAAAAAABcAAAAAAAAACwAAAAAAAAAMAAAAAAAAABUAAAAAAAAAFgAAAAAAAAAMAAAAAAAAAA0AAAAAAAAAFAAAAAAAAAAVAAAAAAAAAA0AAAAAAAAADgAAAAAAAAATAAAAAAAAABQAAAAAAAAADgAAAAAAAAAPAAAAAAAAABIAAAAAAAAAEwAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABAAAAAAAAAAAQAAAAAAAAACAAAAAAAAABEAAAAAAAAAcAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAADgAAAAAAAAAJCQkJCQkJCQkJCQkJCQ==
+  </AppendedData>
+</VTKFile>
diff --git a/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/HC_stabilizier.png b/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/HC_stabilizier.png
new file mode 100644
index 0000000000000000000000000000000000000000..95ad77a37687477e2dc784d2dd5a4cf9ed787f3d
Binary files /dev/null and b/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/HC_stabilizier.png differ
diff --git a/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/index.md b/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/index.md
new file mode 100644
index 0000000000000000000000000000000000000000..c488dffb9f835b94af19797f47763297c40b1b4c
--- /dev/null
+++ b/web/content/docs/benchmarks/hydro-component/Ogata_Banks_isotropic_diffusion_stb/index.md
@@ -0,0 +1,15 @@
++++
+project = ["Parabolic/ComponentTransport/ClassicalTransportExample/classical_transport_example.prj"]
+author = "Wenqing Wang"
+date = "2022-09-12T15:39:05+02:00"
+title = "Ogata Banks' example using the isotropic diffusion stabilization"
+weight = 161
+image = "HC_stabilizier.png"
++++
+
+{{< data-link >}}
+
+See the detailed introduction in
+ [Classical transport example: using the isotropic diffusion stabilization]({{< relref "ClassicalTransportExample_isotropic_diffusion_stb" >}}).
+
+{{< img src="HC_stabilizier.png" >}}