diff --git a/examples/Advection.cpp b/examples/Advection.cpp
new file mode 100644
index 0000000..2ce374f
--- /dev/null
+++ b/examples/Advection.cpp
@@ -0,0 +1,59 @@
+#include "tug/Boundary.hpp"
+#include <cstddef>
+#include <iostream>
+#include <tug/tug.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[]) {
+  constexpr std::size_t row = 5;
+  constexpr std::size_t col = 5;
+
+  RowMajMat<double> hydHeads = RowMajMat<double>::Constant(row, col, 1);
+  RowMajMat<double> concentrations = RowMajMat<double>::Constant(row, col, 0);
+
+  Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+             tug::HYDRAULIC_RESOLVE::EXPLICIT>
+      velocities(hydHeads);
+
+  velocities.setDomain(5, 5);
+  velocities.setPermKX(RowMajMat<double>::Constant(row, col, 3E-7));
+  velocities.setPermKY(RowMajMat<double>::Constant(row, col, 3E-7));
+  velocities.setEpsilon(1E-8);
+
+  Advection advection = Advection(concentrations, velocities);
+
+  advection.setPorosity(RowMajMat<double>::Constant(row, col, 0.2));
+  advection.setIterations(3);
+  // 1 hour
+  advection.setTimestep(1.6666e+06);
+
+  // create boundaries
+  Boundary<double> &bcH = velocities.getBoundaryConditions();
+  bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10);
+  bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcH.setBoundarySideConstant(BC_SIDE_TOP, 1);
+  // bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1);
+  // bcH.setInnerBoundary(row / 2, col / 2, 10);
+
+  Boundary<double> &bcC = advection.getBoundaryConditions();
+  bcC.setBoundarySideConstant(BC_SIDE_LEFT, 1);
+  bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcC.setInnerBoundary(row / 2, col / 2, 1);
+  // bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_TOP, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
+
+  advection.run();
+
+  std::cout << velocities.getConcentrationMatrix() << std::endl << std::endl;
+
+  std::cout << velocities.getVelocitiesX() << std::endl
+            << std::endl
+            << velocities.getVelocitiesY() << std::endl
+            << std::endl;
+
+  std::cout << concentrations << std::endl;
+}
diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp
index c09d93d..f793e5c 100644
--- a/examples/BTCS_2D_proto_example.cpp
+++ b/examples/BTCS_2D_proto_example.cpp
@@ -1,5 +1,5 @@
 #include <Eigen/Eigen>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 3c8bf28..c6f2eb8 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,7 +1,9 @@
-add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
-add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
-add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+# add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
+# add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
+# add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+add_executable(Advection Advection.cpp)
 
-target_link_libraries(BTCS_2D_proto_example tug)
-target_link_libraries(FTCS_2D_proto_closed_mdl tug)
-target_link_libraries(FTCS_2D_proto_example_mdl tug)
+# target_link_libraries(BTCS_2D_proto_example tug)
+# target_link_libraries(FTCS_2D_proto_closed_mdl tug)
+# target_link_libraries(FTCS_2D_proto_example_mdl tug)
+target_link_libraries(Advection tug)
diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp
index 69b31c3..eab6dc9 100644
--- a/examples/FTCS_2D_proto_closed_mdl.cpp
+++ b/examples/FTCS_2D_proto_closed_mdl.cpp
@@ -9,7 +9,7 @@
 #include <Eigen/Eigen>
 #include <cstdlib>
 #include <iostream>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;
diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp
index 41ddb91..f7059a1 100644
--- a/examples/FTCS_2D_proto_example_mdl.cpp
+++ b/examples/FTCS_2D_proto_example_mdl.cpp
@@ -7,7 +7,7 @@
  */
 
 #include <Eigen/Eigen>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;
diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp
new file mode 100644
index 0000000..58ed9b1
--- /dev/null
+++ b/include/tug/Advection/Advection.hpp
@@ -0,0 +1,408 @@
+/**
+ * @file Advection.hpp
+ * @brief API of Advection class, holding information for a simulation of
+ * advection. Holds a predifined Grid object, Boundary object and Velocities
+ * object
+ */
+
+#pragma once
+
+#include "tug/Core/Matrix.hpp"
+#include <Eigen/src/Core/Array.h>
+#include <Eigen/src/Core/AssignEvaluator.h>
+#include <Eigen/src/Core/Matrix.h>
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+#include <tug/Boundary.hpp>
+
+#include <Eigen/Eigen>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
+
+#include <tug/Advection/Velocities.hpp>
+
+using namespace Eigen;
+namespace tug {
+template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
+class Advection : public BaseSimulationGrid<T> {
+private:
+  T timestep{-1};
+  int numThreads{omp_get_num_procs()};
+
+  Velocities<T, hyd_mode, hyd_resolve> &velocities;
+  RowMajMat<T> porosity;
+
+public:
+  /**
+   * @brief Construct a new Advection object, used to calculate material
+   * transport. A timestep and number of iterations must be set. A transient
+   * case can be selected by initializing Steady=false. With each timestep the
+   * Velocities object will also be updated.
+   * A steady case can be selected by initializing Steady=true. The velocities
+   * object will not be updated. Velocities can be simulated to convergence
+   * beforehand. Porosity can be set, the default is 1. CSV Output is off by
+   * default.
+   *
+   * @param grid Valid grid object
+   * @param bc Valid Boundary object
+   * @param Steady Used to choose between Steady and Transient case. Either true
+   * or false
+   */
+  Advection(RowMajMat<T> &origin,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(origin), velocities(_velocities) {
+    tug_assert(origin.rows() == velocities.getConcentrationMatrix().rows() &&
+                   origin.cols() == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
+
+  Advection(T *data, std::size_t rows, std::size_t cols,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(data, rows, cols), velocities(_velocities) {
+    tug_assert(rows == velocities.getConcentrationMatrix().rows() &&
+                   cols == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
+
+  /**
+   * @brief Sets the porosity of the medium
+   *
+   * @param porosity new porosity value
+   */
+  void setPorosity(const RowMajMat<T> &porosity) {
+    tug_assert(porosity.rows() == this->rows() &&
+                   porosity.cols() == this->cols(),
+               "Porosity matrix must have the same dimensions as the grid");
+    tug_assert(porosity.minCoeff() >= 0 && porosity.maxCoeff() <= 1,
+               "Porosity must be a value between 0 and 1 (inclusive)");
+
+    this->porosity = porosity;
+  }
+
+  /**
+   * @brief Set the size of the timestep. Must be greater than zero
+   *
+   * @param timestep Size of the timestep
+   */
+  void setTimestep(T timestep) {
+    tug_assert(timestep > 0, "Timestep must be greater than zero");
+    this->timestep = timestep;
+  }
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Advection
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      if (num_threads > maxThreadNumber) {
+        throw std::invalid_argument(
+            "Number of threads exceeds the number of processor cores (" +
+            std::to_string(maxThreadNumber) + ").");
+      } else {
+        throw std::invalid_argument("Number of threads is less than 1.");
+      }
+    }
+  }
+
+  void run() {
+    this->setDomain(velocities.domainX(), velocities.domainY());
+
+    this->applyInnerBoundaries();
+
+    if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
+      velocities.run();
+    }
+
+    for (int i = 0; i < this->getIterations(); i++) {
+      if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) {
+        velocities.run();
+      }
+      startAdvection();
+    }
+  }
+
+private:
+  /**
+   * @brief function calculating material transport for one timestep
+   */
+  void startAdvection() {
+    const std::size_t rows = this->rows();
+    const std::size_t cols = this->cols();
+    const T dx = this->deltaCol();
+    const T dy = this->deltaRow();
+    const T volume = this->deltaCol() * this->deltaRow();
+
+    RowMajMatMap<T> &newConcentrations = this->getConcentrationMatrix();
+
+    RowMajMat<T> veloX = this->velocities.getVelocitiesX();
+    RowMajMat<T> veloY = this->velocities.getVelocitiesY();
+    const Boundary<T> &bc = this->getBoundaryConditions();
+
+    // Calculate Courant-Levy-Frederich condition
+    const T maxFx = std::max(abs(veloX.maxCoeff()), abs(veloX.minCoeff()));
+    const T maxFy = std::max(abs(veloY.maxCoeff()), abs(veloY.minCoeff()));
+    const T maxFlux = std::max(maxFx, maxFy);
+
+    tug_assert(maxFlux != 0, "Division by zero: maxF is zero.");
+
+    // TODO: sustitute 1 by porisity
+    const RowMajMat<T> porevolume =
+        RowMajMat<T>::Constant(rows, cols, volume) * 1;
+
+    const T cfl = (porevolume / maxFlux).minCoeff();
+    const int innerSteps = (int)ceil(timestep / cfl);
+    const T innerTimestep = timestep / innerSteps;
+
+    std::cout << "CFL (adv) timestep: " << cfl << std::endl;
+    std::cout << "Inner iterations (adv): " << innerSteps << std::endl;
+    std::cout << "Inner timestep (adv): " << innerTimestep << std::endl;
+
+    // const RowMajMat<T> multiplier = porevolume * (1 / innerTimestep);
+
+    RowMajMat<T> fluxX(rows, cols + 1);
+    RowMajMat<T> fluxY(rows + 1, cols);
+
+    for (std::size_t k = 0; k < innerSteps; k++) {
+      const RowMajMat<T> oldConcentrations = newConcentrations;
+
+      // Update flux with concentrations in x-direction
+      for (std::size_t row_i = 0; row_i < fluxX.rows(); row_i++) {
+        for (std::size_t col_i = 0; col_i < fluxX.cols(); col_i++) {
+          T &currentFlux = fluxX(row_i, col_i);
+          const T &currentVelo = veloX(row_i, col_i);
+          const std::int32_t cellIndex = (currentVelo > 0) ? col_i - 1 : col_i;
+
+          if (cellIndex < 0 || cellIndex >= cols) {
+            const auto bcElement = bc.getBoundaryElement(
+                cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i);
+            // if we have a constant boundary, we use the value of the boundary
+            if (bcElement.getType() == BC_TYPE_CONSTANT) {
+              currentFlux = currentVelo * bcElement.getValue();
+              continue;
+            }
+
+            // otherwise it is a closed boundary
+            currentFlux = 0;
+            continue;
+          }
+          // otherwise we use the concentration of the inflow/outflow cell,
+          // multiplied by the velocity
+          currentFlux = currentVelo * oldConcentrations(row_i, cellIndex);
+        }
+      }
+
+      // Update flux with concentrations in y-direction
+      for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) {
+        for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) {
+          T &currentFlux = fluxY(row_i, col_i);
+          const T &currentVelo = veloY(row_i, col_i);
+          const std::int32_t cellIndex = (currentVelo > 0) ? row_i - 1 : row_i;
+
+          if (cellIndex < 0 || cellIndex >= rows) {
+            const auto bcElement = bc.getBoundaryElement(
+                cellIndex < 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i);
+            // if we have a constant boundary, we use the value of the
+            // boundary
+            if (bcElement.getType() == BC_TYPE_CONSTANT) {
+              currentFlux = currentVelo * bcElement.getValue();
+              continue;
+            }
+
+            // otherwise it is a closed boundary
+            currentFlux = 0;
+            continue;
+          }
+          // otherwise we use the concentration of the inflow/outflow cell,
+          // multiplied by the velocity
+          currentFlux = currentVelo * oldConcentrations(cellIndex, col_i);
+        }
+      }
+
+      // Update concentrations
+      for (std::size_t row_i = 0; row_i < rows; row_i++) {
+        for (std::size_t col_i = 0; col_i < cols; col_i++) {
+          const T horizontalFlux =
+              fluxX(row_i, col_i) - fluxX(row_i, col_i + 1);
+          const T verticalFlux = fluxY(row_i, col_i) - fluxY(row_i + 1, col_i);
+          newConcentrations(row_i, col_i) = oldConcentrations(row_i, col_i) +
+                                            innerTimestep /
+                                                porevolume(row_i, col_i) *
+                                                (horizontalFlux + verticalFlux);
+        }
+      }
+
+      // const Boundary<T> leftCorner = bc.getBoundaryElement(BC_SIDE_LEFT, 0);
+      // newConcentrations(0, 0) = oldConcentrations(0, 0) + ;
+      //
+      // for (std::size_t col_i = 0; col_i < cols; col_i++) {
+      //   T horizontalFlux = oldConcentrations(0, col_i);
+      //   if (col_i == 0) {
+      //     switch (bc.getBoundaryElementType(BC_SIDE_LEFT, 0)) {
+      //     case BC_TYPE_CLOSED: {
+      //       break;
+      //     }
+      //     case BC_TYPE_CONSTANT: {
+      //       newConcentrations(0, 0) =
+      //           oldConcentrations(0, 0) +
+      //           (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) * veloX(0, 0)) *
+      //               (2 * innerTimestep) / dx;
+      //       newConcentrations(rows - 1, 0) =
+      //           oldConcentrations(rows - 1, 0) +
+      //           (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) *
+      //            veloX(rows - 1, 0)) *
+      //               (2 * innerTimestep) / dx;
+      //       break;
+      //     }
+      //     }
+      //     newConcentrations(0, 0) -=
+      //         veloX(0, 1) * oldConcentrations(0, 1) * innerTimestep / dx;
+      //   }
+      // }
+      //
+      // for (std::size_t row_i = 1; row_i < rows - 1; row_i++) {
+      //   for (std::size_t col_i = 1; col_i < cols - 1; col_i++) {
+      //     const T fluxHorizontal =
+      //         veloX(row_i, col_i) * oldConcentrations(row_i, col_i - 1) -
+      //         veloX(row_i, col_i + 1) * oldConcentrations(row_i, col_i);
+      //     const T fluxVertical =
+      //         veloY(row_i, col_i) * oldConcentrations(row_i - 1, col_i) -
+      //         veloY(row_i + 1, col_i) * oldConcentrations(row_i, col_i);
+      //     newConcentrations(row_i, col_i) =
+      //         oldConcentrations(row_i, col_i) +
+      //         (fluxHorizontal / dx + fluxVertical / dy) * innerTimestep;
+      //   }
+      // }
+
+      // Calculate sum of incoming/outgoing Flow*Concentration in x-direction in
+      // each cell #pragma omp parallel for num_threads(numThreads)
+      //       for (int i = 0; i < rows; i++) {
+      //         for (int j = 0; j < cols + 1; j++) {
+      //           if (j == 0) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloX(i, j) > 0) {
+      //                 // outx positive -> flow from border to cell i,j
+      //                 newConcentrations(i, j) +=
+      //                     veloX(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+      //               } else if (veloX(i, j) < 0) {
+      //                 // outx negative -> flow from i,j towards border
+      //                 newConcentrations(i, j) += veloX(i, j) *
+      //                 oldConcentrations(i, j);
+      //               }
+      //             }
+      //           } else if (j == cols) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloX(i, j) > 0) {
+      //                 // outx positive-> flow from i,j-1 towards border
+      //                 newConcentrations(i, j - 1) -=
+      //                     veloX(i, j) * oldConcentrations(i, j - 1);
+      //               } else if (veloX(i, j) < 0) {
+      //                 // outx negative -> flow from border to cell i,j-1
+      //                 newConcentrations(i, j - 1) -=
+      //                     veloX(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+      //               }
+      //             }
+      //           }
+      //           // flow between inner cells
+      //           else {
+      //             // outx positive -> flow from cell i,j-1 towards cell i,j
+      //             if (veloX(i, j) > 0) {
+      //               newConcentrations(i, j - 1) -=
+      //                   veloX(i, j) * oldConcentrations(i, j - 1);
+      //               newConcentrations(i, j) +=
+      //                   veloX(i, j) * oldConcentrations(i, j - 1);
+      //             }
+      //             // outx negative -> flow from cell i,j toward cell i,j-1
+      //             else if (veloX(i, j) < 0) {
+      //               newConcentrations(i, j - 1) +=
+      //                   veloX(i, j) * oldConcentrations(i, j);
+      //               newConcentrations(i, j) += veloX(i, j) *
+      //               oldConcentrations(i, j);
+      //             }
+      //           }
+      //         }
+      //       }
+      // // calculate sum in y-direction
+      // // parallelize outer loop over columns to ensure thread-safety, each
+      // thread only
+      // // modifies cells within its column
+      // #pragma omp parallel for num_threads(numThreads)
+      //       for (int j = 0; j < cols; j++) {
+      //         for (int i = 0; i < rows + 1; i++) {
+      //           if (i == 0) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_TOP, j) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloY(i, j) > 0) {
+      //                 // outy positive -> flow from border to cell i,j
+      //                 newConcentrations(i, j) +=
+      //                     veloY(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_TOP, j);
+      //               } else if (veloY(i, j) < 0) {
+      //                 // outy negative -> flow from i,j towards border
+      //                 newConcentrations(i, j) += veloY(i, j) *
+      //                 oldConcentrations(i, j);
+      //               }
+      //             }
+      //           } else if (i == rows) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) !=
+      //                 BC_TYPE_CLOSED) {
+      //               if (veloY(i, j) > 0) {
+      //                 // outy positive-> flow from i-1,j towards border
+      //                 newConcentrations(i - 1, j) -=
+      //                     veloY(i, j) * oldConcentrations(i - 1, j);
+      //               } else if (veloY(i, j) < 0) {
+      //                 // outy negative -> flow from border to cell i,j-1
+      //                 newConcentrations(i - 1, j) -=
+      //                     veloY(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j);
+      //               }
+      //             }
+      //           }
+      //           // flow between inner cells
+      //           else {
+      //             // outy positive -> flow from cell i-1,j towards cell i,j
+      //             if (veloY(i, j) > 0) {
+      //               newConcentrations(i - 1, j) -=
+      //                   veloY(i, j) * oldConcentrations(i - 1, j);
+      //               newConcentrations(i, j) +=
+      //                   veloY(i, j) * oldConcentrations(i - 1, j);
+      //             }
+      //             // outy negative -> flow from cell i,j toward cell i-1,j
+      //             else if (veloY(i, j) < 0) {
+      //               newConcentrations(i - 1, j) -=
+      //                   veloY(i, j) * oldConcentrations(i, j);
+      //               newConcentrations(i, j) += veloY(i, j) *
+      //               oldConcentrations(i, j);
+      //             }
+      //           }
+      //         }
+      //       }
+
+      // for (std::size_t row_i = 0; row_i < rows; row_i++) {
+      //   for (std::size_t col_i = 0; col_i < cols; col_i++) {
+      //     newConcentrations(row_i, col_i) =
+      //         oldConcentrations(row_i, col_i) +
+      //         newConcentrations(row_i, col_i) * multiplier(row_i, col_i);
+      //   }
+      // }
+    }
+  }
+};
+} // namespace tug
diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp
new file mode 100644
index 0000000..6679db8
--- /dev/null
+++ b/include/tug/Advection/Velocities.hpp
@@ -0,0 +1,376 @@
+/**
+ * @file Velocities.hpp
+ * @brief API of Velocities class, holding information for a simulation of
+ * Hydraulic Charge and Darcy-Velocities. Holds a predifined Grid object and
+ * Boundary object.
+ *
+ */
+
+#pragma once
+
+#include "tug/Core/Numeric/SimulationInput.hpp"
+#include <algorithm>
+#include <cstddef>
+#include <iostream>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+
+#include <tug/Boundary.hpp>
+#include <tug/Core/BaseSimulation.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_num_procs() 1
+#endif
+
+using namespace Eigen;
+namespace tug {
+
+enum class HYDRAULIC_MODE { TRANSIENT, STEADY_STATE };
+enum class HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT };
+
+template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
+class Velocities : public BaseSimulationGrid<T> {
+private:
+  int innerIterations{1};
+  T timestep{-1};
+  T epsilon{1E-5};
+  int numThreads{omp_get_num_procs()};
+  bool steady{false};
+
+  RowMajMat<T> velocitiesX;
+  RowMajMat<T> velocitiesY;
+
+  RowMajMat<T> permKX;
+  RowMajMat<T> permKY;
+
+public:
+  Velocities(RowMajMat<T> &origin)
+      : BaseSimulationGrid<T>(origin),
+        velocitiesX(origin.rows(), origin.cols() + 1),
+        velocitiesY(origin.rows() + 1, origin.cols()),
+        permKX(origin.rows(), origin.cols()),
+        permKY(origin.rows(), origin.cols()) {};
+
+  Velocities(T *data, std::size_t rows, std::size_t cols)
+      : BaseSimulationGrid<T>(data, rows, cols), velocitiesX(rows, cols + 1),
+        velocitiesY(rows + 1, cols), permKX(rows, cols), permKY(rows, cols) {};
+
+  // Velocities(T *data, std::size_t length)
+  //     : BaseSimulationGrid<T>(data, 1, length), velocitiesX(1, length + 1),
+  //       alphaX(1, length) {};
+
+  /**
+   * @brief Set the epsilon value, the relativ error allowed for convergence
+   *
+   * @param epsilon the new epsilon value
+   */
+  void setEpsilon(T epsilon) {
+    tug_assert(0 <= epsilon && epsilon < 1,
+               "Relative Error epsilon must be between 0 and 1");
+
+    this->epsilon = epsilon;
+  }
+
+  /**
+   * @brief Get the alphaX matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaX matrix.
+   */
+  RowMajMat<T> &getPermKX() { return permKX; }
+
+  /**
+   * @brief Get the alphaY matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaY matrix.
+   */
+  RowMajMat<T> &getPermKY() {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    return permKY;
+  }
+
+  /**
+   * @brief Set the alphaX matrix.
+   *
+   * @param alphaX The new alphaX matrix.
+   */
+  void setPermKX(const RowMajMat<T> &alphaX) { this->permKX = alphaX; }
+
+  /**
+   * @brief Set the alphaY matrix.
+   *
+   * @param alphaY The new alphaY matrix.
+   */
+  void setPermKY(const RowMajMat<T> &alphaY) {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    this->permKY = alphaY;
+  }
+
+  bool isSteady() const { return steady; }
+
+  /**
+   * @brief Set the timestep per iteration
+   *
+   * @param timestep timestep per iteration
+   */
+  void setTimestep(T timestep) override {
+    if (timestep <= 0) {
+      throw std::invalid_argument("Timestep must be greater than zero");
+    }
+    this->timestep = timestep;
+    const T deltaColSquare = this->deltaCol() * this->deltaCol();
+    const T deltaRowSquare = this->deltaRow() * this->deltaRow();
+    const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+
+    const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff());
+
+    T cfl = minDeltaSquare / (4 * maxK);
+    if (timestep > cfl) {
+      this->innerIterations = (int)ceil(timestep / cfl);
+      this->timestep = timestep / (double)innerIterations;
+      std::cerr << "Warning :: Timestep was adjusted, because of stability "
+                   "conditions. Time duration was approximately preserved by "
+                   "adjusting internal number of iterations."
+                << std::endl;
+      std::cout << "FTCS" << "_" << "2D" << " :: Required "
+                << this->innerIterations
+                << " inner iterations with dt=" << this->timestep << std::endl;
+    } else {
+      this->innerIterations = 1;
+    }
+  };
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Velocities
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      throw std::invalid_argument(
+          "Number of threads exceeds the number of processor cores (" +
+          std::to_string(maxThreadNumber) + ") or is less than 1.");
+    }
+  }
+
+  /**
+   * @brief Getter function for outx, the matrix containing velocities in
+   * x-Direction; returns a reference to outx
+   *
+   * */
+  const RowMajMat<T> &getVelocitiesX() const { return this->velocitiesX; }
+
+  /**
+   * @brief Getter function for outy, the matrix containing velocities in
+   * y-Direction; return a reference to outy
+   */
+  const RowMajMat<T> &getVelocitiesY() const { return this->velocitiesY; }
+
+  /**
+   * @brief Simulation of hydraulic charge either until convergence,
+   * or for a number of selected timesteps. Calculation of Darcy-velocities.
+   */
+  void run() override {
+    // if iterations < 1 calculate hydraulic charge until steady state is
+    // reached
+
+    this->applyInnerBoundaries();
+
+    SimulationInput<T> input = {.concentrations =
+                                    this->getConcentrationMatrix(),
+                                .alphaX = this->getPermKX(),
+                                .alphaY = this->getPermKY(),
+                                .boundaries = this->getBoundaryConditions(),
+                                .dim = this->getDim(),
+                                .timestep = this->timestep,
+                                .rowMax = this->rows(),
+                                .colMax = this->cols(),
+                                .deltaRow = this->deltaRow(),
+                                .deltaCol = this->deltaCol()};
+
+    if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
+      if (steady) {
+        return;
+      }
+
+      const T deltaColSquare = this->deltaCol() * this->deltaCol();
+      const T deltaRowSquare = this->deltaRow() * this->deltaRow();
+      const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+
+      const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff());
+      // Calculate largest possible timestep, depending on K and gridsize
+
+      std::cout << "CFL (hydHead) timestep: " << minDeltaSquare / (4 * maxK)
+                << std::endl;
+
+      setTimestep(minDeltaSquare / (4 * maxK));
+
+      input.timestep = this->timestep;
+
+      RowMajMat<T> oldConcentrations;
+      do {
+        oldConcentrations = this->getConcentrationMatrix();
+        (void)calculate_hydraulic_flow(input);
+      } while (!checkConvergance(oldConcentrations));
+
+      steady = true;
+
+    } else {
+      if (timestep == -1) {
+        throw_invalid_argument("Timestep is not set");
+      }
+      for (int i = 0; i < innerIterations; i++) {
+        (void)calculate_hydraulic_flow(input);
+      }
+    }
+
+    (void)computeFluidVelocities();
+  };
+
+private:
+  /**
+   * @brief Calculate the new hydraulic charge using FTCS
+   */
+  void calculate_hydraulic_flow(SimulationInput<T> &sim_in) {
+    if constexpr (hyd_resolve == HYDRAULIC_RESOLVE::EXPLICIT) {
+      FTCS_2D(sim_in, numThreads);
+    } else {
+      BTCS_2D(sim_in, ThomasAlgorithm, numThreads);
+    }
+  };
+
+  /**
+   * @brief checks if the matrix of Hydraulic Heads has converged
+   *
+   * @return bool true if for all corresponding cells of the matrices,
+   * containing old and new Charge values, the relative error is below the
+   * selected Epsilon
+   */
+  bool checkConvergance(const RowMajMat<T> &oldHeads) {
+    const auto abs_err = (oldHeads - this->getConcentrationMatrix()).cwiseAbs();
+    const auto rel_err = abs_err.cwiseQuotient(this->getConcentrationMatrix());
+
+    return rel_err.maxCoeff() < epsilon;
+  }
+
+  /**
+   * @brief Update the matrices containing Darcy velocities in x and y
+   * directions
+   */
+  void computeFluidVelocities() {
+    const std::size_t rows = this->rows();
+    const std::size_t cols = this->cols();
+    const T dx = this->deltaCol();
+    const T dy = this->deltaRow();
+    const RowMajMat<T> &hydraulicCharges = this->getConcentrationMatrix();
+
+    const RowMajMat<T> &permKX = this->permKX;
+    const RowMajMat<T> &permKY = this->permKY;
+
+    const Boundary<T> &bc = this->getBoundaryConditions();
+
+    // calculate velocities in x-direction
+    for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+      const auto bc_left = bc.getBoundaryElement(BC_SIDE_LEFT, i_rows);
+      switch (bc_left.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesX(i_rows, 0) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesX(i_rows, 0) =
+            -this->permKX(i_rows, 0) *
+            (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2);
+        break;
+      }
+      }
+
+      const auto bc_right = bc.getBoundaryElement(BC_SIDE_RIGHT, i_rows);
+      switch (bc_right.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesX(i_rows, cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesX(i_rows, cols) =
+            -permKX(i_rows, cols - 1) *
+            (bc_right.getValue() - hydraulicCharges(i_rows, cols - 1)) /
+            (dx / 2);
+        break;
+      }
+      }
+    }
+
+// main loop for calculating velocities in x-direction for inner cells
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 0; i < rows; i++) {
+      for (int j = 1; j < cols; j++) {
+        velocitiesX(i, j) =
+            -permKX(i, j - 1) *
+            (hydraulicCharges(i, j) - hydraulicCharges(i, j - 1)) / dx;
+      }
+    }
+
+    // calculate velocities in y-direction
+    for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+      const auto bc_top = bc.getBoundaryElement(BC_SIDE_TOP, i_cols);
+      switch (bc_top.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesY(0, i_cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesY(0, i_cols) =
+            -permKY(0, i_cols) *
+            (hydraulicCharges(0, i_cols) - bc_top.getValue()) / (dy / 2);
+        break;
+      }
+      }
+
+      const auto bc_bottom = bc.getBoundaryElement(BC_SIDE_BOTTOM, i_cols);
+      switch (bc_bottom.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesY(rows, i_cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesY(rows, i_cols) =
+            -permKY(rows - 1, i_cols) *
+            (bc_bottom.getValue() - hydraulicCharges(rows - 1, i_cols)) /
+            (dy / 2);
+        break;
+      }
+      }
+    }
+
+// main loop for calculating velocities in y-direction for inner cells
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 1; i < rows; i++) {
+      for (int j = 0; j < cols; j++) {
+        velocitiesY(i, j) =
+            -permKY(i, j) *
+            (hydraulicCharges(i, j) - hydraulicCharges(i - 1, j)) / dy;
+      }
+    }
+  };
+};
+} // namespace tug
diff --git a/include/tug/Core/BaseSimulation.hpp b/include/tug/Core/BaseSimulation.hpp
index 82a1480..3dc8268 100644
--- a/include/tug/Core/BaseSimulation.hpp
+++ b/include/tug/Core/BaseSimulation.hpp
@@ -64,6 +64,18 @@ private:
   T delta_col;
   T delta_row;
 
+protected:
+  void applyInnerBoundaries() {
+    const auto &inner_bc = boundaryConditions.getInnerBoundaries();
+    if (inner_bc.empty()) {
+      return;
+    }
+
+    for (const auto &[rowcol, value] : inner_bc) {
+      concentrationMatrix(rowcol.first, rowcol.second) = value;
+    }
+  }
+
 public:
   /**
    * @brief Constructs a BaseSimulationGrid from a given RowMajMat object.
diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp
index ce1d5e8..4ccf1a4 100644
--- a/include/tug/Core/Numeric/FTCS.hpp
+++ b/include/tug/Core/Numeric/FTCS.hpp
@@ -8,12 +8,13 @@
 #ifndef FTCS_H_
 #define FTCS_H_
 
-#include "tug/Core/TugUtils.hpp"
+#include "tug/Core/Matrix.hpp"
 #include <cstddef>
 #include <cstring>
 #include <tug/Boundary.hpp>
 #include <tug/Core/Matrix.hpp>
 #include <tug/Core/Numeric/SimulationInput.hpp>
+#include <tug/Core/TugUtils.hpp>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -53,6 +54,7 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
   }
 
   tug_assert(false, "Undefined Boundary Condition Type!");
+  return 0;
 }
 
 // FTCS solution for 1D grid
@@ -62,18 +64,23 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   const T &timestep = input.timestep;
 
   RowMajMatMap<T> &concentrations_grid = input.concentrations;
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
   const auto &alphaX = input.alphaX;
   const auto &bc = input.boundaries;
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
 
+  const auto inner_bc = bc.getInnerBoundaryRow(0);
+
   // inner cells
   // independent of boundary condition type
   for (int col = 1; col < colMax - 1; col++) {
+    if (inner_bc[col].first) {
+      continue;
+    }
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_left = concentrations_grid(row, col - 1);
     const T &conc_right = concentrations_grid(row, col + 1);
@@ -90,8 +97,8 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   }
 
   // left boundary; hold column constant at index 0
-  {
-    int col = 0;
+  int col = 0;
+  if (!inner_bc[col].first) {
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_right = concentrations_grid(row, col + 1);
     const T &alpha_c = alphaX(row, col);
@@ -107,8 +114,8 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   }
 
   // right boundary; hold column constant at max index
-  {
-    int col = colMax - 1;
+  col = colMax - 1;
+  if (!inner_bc[col].first) {
     const T &conc_c = concentrations_grid(row, col);
     const T &conc_left = concentrations_grid(row, col - 1);
     const T &alpha_c = alphaX(row, col);
@@ -136,10 +143,6 @@ static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
   const T &timestep = input.timestep;
 
   RowMajMatMap<T> &concentrations_grid = input.concentrations;
-
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = concentrations_grid;
-
   const auto &alphaX = input.alphaX;
   const auto &alphaY = input.alphaY;
   const auto &bc = input.boundaries;
@@ -147,13 +150,18 @@ static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
   const T sx = timestep / (deltaCol * deltaCol);
   const T sy = timestep / (deltaRow * deltaRow);
 
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
+
 #pragma omp parallel for num_threads(numThreads)
   for (std::size_t row_i = 0; row_i < rowMax; row_i++) {
     for (std::size_t col_i = 0; col_i < colMax; col_i++) {
+      if (bc.getInnerBoundary(row_i, col_i).first) {
+        continue;
+      }
       // horizontal change
       T horizontal_change;
       {
-
         const T &conc_c = concentrations_grid(row_i, col_i);
         const T &alpha_c = alphaX(row_i, col_i);
 
diff --git a/include/tug/Core/Numeric/SimulationInput.hpp b/include/tug/Core/Numeric/SimulationInput.hpp
index 202619c..5068537 100644
--- a/include/tug/Core/Numeric/SimulationInput.hpp
+++ b/include/tug/Core/Numeric/SimulationInput.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <cstdint>
 #include <tug/Boundary.hpp>
 #include <tug/Core/Matrix.hpp>
 
@@ -10,12 +11,11 @@ template <typename T> struct SimulationInput {
   const RowMajMat<T> &alphaX;
   const RowMajMat<T> &alphaY;
   const Boundary<T> boundaries;
-
   const std::uint8_t dim;
-  const T timestep;
+  T timestep;
   const std::size_t rowMax;
   const std::size_t colMax;
   const T deltaRow;
   const T deltaCol;
 };
-} // namespace tug
\ No newline at end of file
+} // namespace tug
diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion/Diffusion.hpp
similarity index 99%
rename from include/tug/Diffusion.hpp
rename to include/tug/Diffusion/Diffusion.hpp
index e353550..f977514 100644
--- a/include/tug/Diffusion.hpp
+++ b/include/tug/Diffusion/Diffusion.hpp
@@ -347,6 +347,8 @@ public:
 
     auto begin = std::chrono::high_resolution_clock::now();
 
+    this->applyInnerBoundaries();
+
     SimulationInput<T> sim_input = {.concentrations =
                                         this->getConcentrationMatrix(),
                                     .alphaX = this->getAlphaX(),
diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp
new file mode 100644
index 0000000..6770de3
--- /dev/null
+++ b/include/tug/tug.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+#include <tug/Advection/Advection.hpp>
+#include <tug/Advection/Velocities.hpp>
+#include <tug/Boundary.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 2a8d780..c0d2a51 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -12,6 +12,8 @@ FetchContent_MakeAvailable(googletest)
 add_executable(testTug 
 setup.cpp
 testDiffusion.cpp 
+testVelocities.cpp
+testAdvection.cpp
 testFTCS.cpp 
 testBoundary.cpp
 )
diff --git a/test/testAdvection.cpp b/test/testAdvection.cpp
new file mode 100644
index 0000000..ac4c2db
--- /dev/null
+++ b/test/testAdvection.cpp
@@ -0,0 +1,69 @@
+#include <gtest/gtest.h>
+#include <tug/tug.hpp>
+
+#define ADVECTION_TEST(x) TEST(Advection, x)
+
+ADVECTION_TEST(LeftToRight) {
+  constexpr std::size_t rows = 21;
+  constexpr std::size_t cols = 21;
+
+  constexpr double K = 1E-2;
+  constexpr double timestep = 5039.05;
+  constexpr std::size_t iterations = 21;
+  constexpr double porosity = 0.2;
+
+  constexpr double epsilon = 1E-13;
+
+  tug::RowMajMat<double> hydHeads =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+  tug::RowMajMat<double> concentrations =
+      tug::RowMajMat<double>::Constant(rows, cols, 0);
+
+  tug::RowMajMat<double> permK =
+      tug::RowMajMat<double>::Constant(rows, cols, K);
+
+  tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+                  tug::HYDRAULIC_RESOLVE::IMPLICIT>
+      velocities(hydHeads);
+  velocities.setDomain(100, 100);
+  velocities.setPermKX(permK);
+  velocities.setPermKY(permK);
+  velocities.setEpsilon(1E-8);
+
+  tug::Advection advection(concentrations, velocities);
+
+  advection.setPorosity(tug::RowMajMat<double>::Constant(rows, cols, porosity));
+  advection.setIterations(iterations);
+  advection.setTimestep(timestep);
+
+  tug::Boundary<double> &bcH = velocities.getBoundaryConditions();
+  bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 10);
+  bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
+
+  tug::Boundary<double> &bcC = advection.getBoundaryConditions();
+  bcC.setBoundarySideConstant(tug::BC_SIDE_LEFT, 0.1);
+  bcC.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
+
+  advection.run();
+
+  // check if the concentration is transported from left to right
+  for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols - 1; i_cols++) {
+      if (i_cols == 0) {
+        EXPECT_LE(concentrations(i_rows, i_cols), 10);
+      } else {
+        EXPECT_GE(concentrations(i_rows, i_cols),
+                  concentrations(i_rows, i_cols + 1));
+      }
+    }
+  }
+
+  // the values should also be equal from top to bottom
+  for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+    const double &ref = concentrations(0, i_cols);
+    for (std::size_t i_rows = 1; i_rows < rows; i_rows++) {
+      // check if the values are equal within the epsilon range
+      EXPECT_NEAR(ref, concentrations(i_rows, i_cols), epsilon);
+    }
+  }
+}
diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp
index 1a7c11a..99e8668 100644
--- a/test/testDiffusion.cpp
+++ b/test/testDiffusion.cpp
@@ -3,7 +3,7 @@
 #include "gtest/gtest.h"
 #include <gtest/gtest.h>
 #include <stdexcept>
-#include <tug/Diffusion.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
 
 #include <Eigen/src/Core/Matrix.h>
 #include <string>
@@ -243,3 +243,47 @@ DIFFUSION_TEST(ConstantInnerCell) {
 
   EXPECT_FALSE((concentrations_result.array() < 0.0).any());
 }
+
+DIFFUSION_TEST(Symmetry) {
+  // Arrange
+  constexpr std::size_t rows = 25;
+  constexpr std::size_t cols = 25;
+
+  constexpr std::size_t center_row = rows / 2;
+  constexpr std::size_t center_col = cols / 2;
+
+  tug::RowMajMat<double> concentrations =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+
+  tug::RowMajMat<double> alpha =
+      tug::RowMajMat<double>::Constant(rows, cols, 1E-2);
+
+  tug::Diffusion<double, tug::FTCS_APPROACH> sim(concentrations);
+
+  sim.setDomain(100, 100);
+  sim.setAlphaX(alpha);
+  sim.setAlphaY(alpha);
+  // choose a high number of iterations, which lead to small changes in ULP
+  // between symmetric cells
+  sim.setIterations(10000);
+  sim.setTimestep(10);
+
+  tug::Boundary<double> &bcH = sim.getBoundaryConditions();
+
+  bcH.setInnerBoundary(center_row, center_col, 10);
+
+  sim.run();
+
+  // check symmetry
+  for (std::size_t i_rows = 0; i_rows <= center_row; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols <= center_col; i_cols++) {
+      if (i_rows == center_row && i_cols == center_col) {
+        continue;
+      }
+      // to avoid floating point errors, we check with ASSERT_DOUBLE_EQ with a
+      // precision of ULP(4), see https://stackoverflow.com/a/4149599
+      ASSERT_DOUBLE_EQ(concentrations(i_rows, i_cols),
+                       concentrations(rows - i_rows - 1, cols - i_cols - 1));
+    }
+  }
+}
diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp
new file mode 100644
index 0000000..40f16d2
--- /dev/null
+++ b/test/testVelocities.cpp
@@ -0,0 +1,78 @@
+#include "tug/Boundary.hpp"
+#include "tug/Core/Matrix.hpp"
+#include <cstddef>
+#include <tug/Advection/Velocities.hpp>
+
+#include <gtest/gtest.h>
+
+#define VELOCITIES_TEST(x) TEST(Velocities, x)
+
+VELOCITIES_TEST(SteadyStateCenter) {
+  // Arrange
+  constexpr std::size_t rows = 25;
+  constexpr std::size_t cols = 25;
+
+  constexpr std::size_t center_row = rows / 2;
+  constexpr std::size_t center_col = cols / 2;
+
+  constexpr double K = 1E-2;
+
+  tug::RowMajMat<double> hydHeads =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+
+  tug::RowMajMat<double> permK =
+      tug::RowMajMat<double>::Constant(rows, cols, K);
+
+  tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+                  tug::HYDRAULIC_RESOLVE::EXPLICIT>
+      velo(hydHeads);
+
+  velo.setDomain(100, 100);
+  velo.setPermKX(permK);
+  velo.setPermKY(permK);
+
+  tug::Boundary<double> &bcH = velo.getBoundaryConditions();
+
+  bcH.setInnerBoundary(center_row, center_col, 10);
+
+  velo.run();
+
+  const auto &velocitiesX = velo.getVelocitiesX();
+  const auto &velocitiesY = velo.getVelocitiesY();
+
+  // Assert
+
+  // check velocities in x-direction
+  for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols + 1; i_cols++) {
+      if (i_rows <= center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows <= center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesX(i_rows, i_cols), 0);
+      } else {
+        FAIL() << "Uncovered case";
+      }
+    }
+  }
+
+  // check velocities in y-direction
+  for (std::size_t i_rows = 0; i_rows < rows + 1; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+      if (i_rows <= center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows <= center_row && i_cols > center_col) {
+        EXPECT_LE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols <= center_col) {
+        EXPECT_GE(velocitiesY(i_rows, i_cols), 0);
+      } else {
+        FAIL() << "Uncovered case";
+      }
+    }
+  }
+}