diff --git a/Applications/ApplicationsLib/CMakeLists.txt b/Applications/ApplicationsLib/CMakeLists.txt
index 701c66a4c58918e3921934ed248a42f9ba381ba7..64bb9c8717445ea5e5994728959e5d14fd50af23 100644
--- a/Applications/ApplicationsLib/CMakeLists.txt
+++ b/Applications/ApplicationsLib/CMakeLists.txt
@@ -10,6 +10,10 @@ target_link_libraries(ApplicationsLib
     PRIVATE MathLib MeshLib MeshGeoToolsLib NumLib
 )
 
+if(OGS_USE_PYTHON)
+    target_include_directories(ApplicationsLib PRIVATE ${PYTHON_INCLUDE_DIRS})
+endif()
+
 # Set cpp definitions if the cmake option is enabled for the given process.
 foreach(process ${ProcessesList})
     if(OGS_BUILD_PROCESS_${process})
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 58b077d4a0af2460de993c2f215b890aaed47997..de782c13f956efde4b975147ad9a0efd64d8fa6d 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -18,6 +18,10 @@
 
 #include <logog/include/logog.hpp>
 
+#ifdef OGS_USE_PYTHON
+#include <pybind11/eval.h>
+#endif
+
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/FileTools.h"
 
@@ -162,6 +166,21 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
 {
     _mesh_vec = readMeshes(project_config, project_directory);
 
+#ifdef OGS_USE_PYTHON
+    if (auto const python_script =
+            project_config.getConfigParameterOptional<std::string>(
+                "python_script"))
+    {
+        namespace py = pybind11;
+        auto const script_path =
+            BaseLib::copyPathToFileName(*python_script, project_directory);
+
+        // Evaluate in scope of main module
+        py::object scope = py::module::import("__main__").attr("__dict__");
+        py::eval_file(script_path, scope);
+    }
+#endif  // OGS_USE_PYTHON
+
     //! \ogs_file_param{prj__curves}
     parseCurves(project_config.getConfigSubtreeOptional("curves"));
 
diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index c2410af0547a690a67d35040746b99cf947c35b4..998ac84220bb3144c299669aaaf317c5674b1763 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/Applications/CLI/CMakeLists.txt
@@ -5,6 +5,46 @@ target_link_libraries(ogs
     PRIVATE BaseLib ApplicationsLib
 )
 
+if(OGS_USE_PYTHON)
+    # Troubleshooting:
+    # If you get linker errors, such as   ogs.cpp:(.text+0xb4): Warnung: undefinierter Verweis auf »_Py_ZeroStruct«
+    # it could be that OGS is compiled with the wrong Python version.
+    # I (Ch. Leh.) observed the following: The symbol _Py_ZeroStruct could not be found in /usr/lib/libpython3.6m.so (I intended to compile OGS with Python3).
+    # It's apparently a Python2 symbol (present in /usr/lib/libpython2.7.so)
+    # The compiler command-line was the following:
+    # /usr/bin/g++ ... -DvtkRenderingVolume_AUTOINIT="1(vtkRenderingVolumeOpenGL2)" \
+    #     -I/usr/include/vtk -I/usr/include/python2.7 -I/usr/include/freetype2 \
+    #     -I/usr/include/libxml2 ... -I/.../BaseLib ... \
+    #     -isystem /usr/include/python3.6m ... -o CMakeFiles/ogs.dir/ogs.cpp.o \
+    #     -c /.../Applications/CLI/ogs.cpp
+    # In particular, the Python2 include path comes before the Python3 include path.
+    # Compiling OGS with Python2 solved the issue.
+    # I assume (this is only a guess!) that VTK pulls in Python2 dependencies (on my system).
+    # I assume that this is related to https://github.com/ufz/ogs/pull/2158.
+    # Workaround: Always make sure that OGS is compiled with the same Python version as VTK.
+    # TODO: Find out how to properly address the issue.
+
+    # Performance warning from
+    # https://github.com/pybind/pybind11/blob/master/docs/compiling.rst:
+    # Since pybind11 is a metatemplate library, it is crucial that certain compiler
+    # flags are provided to ensure high quality code generation. In contrast to the
+    # pybind11_add_module() command, the CMake interface library only provides the
+    # minimal set of parameters to ensure that the code using pybind11 compiles, but
+    # it does not pass these extra compiler flags (i.e. this is up to you).
+    # TODO: Enable further compiler/linker flags.
+    target_link_libraries(ogs PRIVATE pybind11::embed)
+
+    add_library(ogs_python_bindings STATIC ogs_python_bindings.cpp)
+    target_link_libraries(ogs_python_bindings PRIVATE pybind11::embed)
+    target_include_directories(ogs_python_bindings PRIVATE ${PYTHON_INCLUDE_DIRS})
+
+    target_link_libraries(ogs PRIVATE ogs_python_bindings)
+
+    if(BUILD_SHARED_LIBS)
+        target_compile_definitions(ogs_python_bindings OGS_BUILD_SHARED_LIBS)
+    endif()
+endif()
+
 if(OGS_USE_PETSC)
     target_link_libraries(ogs PRIVATE ${PETSC_LIBRARIES})
 endif()
diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index 1a7dbc096a267bedfe71c3135d74e77dfddc1e36..f9ceb6a55bac2e23c5598f815f2be3788559b749 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -10,8 +10,8 @@
  *
  */
 
-#include <chrono>
 #include <tclap/CmdLine.h>
+#include <chrono>
 
 #ifndef _WIN32
 #ifdef __APPLE__
@@ -42,7 +42,9 @@
 
 #include "NumLib/NumericsConfig.h"
 
-int main(int argc, char *argv[])
+#include "ogs_python_bindings.h"
+
+int main(int argc, char* argv[])
 {
     // Parse CLI arguments.
     TCLAP::CmdLine cmd(
@@ -65,34 +67,32 @@ int main(int argc, char *argv[])
         "PROJECT FILE");
     cmd.add(project_arg);
 
-    TCLAP::ValueArg<std::string> outdir_arg(
-        "o", "output-directory",
-        "the output directory to write to",
-        false,
-        "",
-        "output directory");
+    TCLAP::ValueArg<std::string> outdir_arg("o", "output-directory",
+                                            "the output directory to write to",
+                                            false, "", "output directory");
     cmd.add(outdir_arg);
 
-    TCLAP::ValueArg<std::string> log_level_arg(
-        "l", "log-level",
-        "the verbosity of logging messages: none, error, warn, info, debug, all",
-        false,
+    TCLAP::ValueArg<std::string> log_level_arg("l", "log-level",
+                                               "the verbosity of logging "
+                                               "messages: none, error, warn, "
+                                               "info, debug, all",
+                                               false,
 #ifdef NDEBUG
-        "info",
+                                               "info",
 #else
-        "all",
+                                               "all",
 #endif
-        "log level");
+                                               "log level");
     cmd.add(log_level_arg);
 
     TCLAP::SwitchArg nonfatal_arg("",
-        "config-warnings-nonfatal",
-        "warnings from parsing the configuration file will not trigger program abortion");
+                                  "config-warnings-nonfatal",
+                                  "warnings from parsing the configuration "
+                                  "file will not trigger program abortion");
     cmd.add(nonfatal_arg);
 
-    TCLAP::SwitchArg unbuffered_cout_arg("",
-        "unbuffered-std-out",
-        "use unbuffered standard output");
+    TCLAP::SwitchArg unbuffered_cout_arg("", "unbuffered-std-out",
+                                         "use unbuffered standard output");
     cmd.add(unbuffered_cout_arg);
 
 #ifndef _WIN32  // TODO: On windows floating point exceptions are not handled
@@ -124,6 +124,11 @@ int main(int argc, char *argv[])
 #endif  // __APPLE__
 #endif  // _WIN32
 
+#ifdef OGS_USE_PYTHON
+    pybind11::scoped_interpreter guard = ApplicationsLib::setupEmbeddedPython();
+    (void)guard;
+#endif
+
     BaseLib::RunTime run_time;
 
     {
@@ -167,7 +172,9 @@ int main(int argc, char *argv[])
             if (auto t = project_config->getConfigSubtreeOptional("insitu"))
             {
                 //! \ogs_file_param{prj__insitu__scripts}
-                InSituLib::Initialize(t->getConfigSubtree("scripts"), BaseLib::extractPath(project_arg.getValue()));
+                InSituLib::Initialize(
+                    t->getConfigSubtree("scripts"),
+                    BaseLib::extractPath(project_arg.getValue()));
                 isInsituConfigured = true;
             }
 #else
@@ -196,11 +203,11 @@ int main(int argc, char *argv[])
 #ifdef USE_INSITU
             if (isInsituConfigured)
                 InSituLib::Finalize();
- #endif
+#endif
             INFO("[time] Execution took %g s.", run_time.elapsed());
 
 #if defined(USE_PETSC)
-            controller->Finalize(1) ;
+            controller->Finalize(1);
 #endif
         }  // This nested scope ensures that everything that could possibly
            // possess a ConfigTree is destructed before the final check below is
@@ -209,7 +216,9 @@ int main(int argc, char *argv[])
         BaseLib::ConfigTree::assertNoSwallowedErrors();
 
         ogs_status = solver_succeeded ? EXIT_SUCCESS : EXIT_FAILURE;
-    } catch (std::exception& e) {
+    }
+    catch (std::exception& e)
+    {
         ERR(e.what());
         ogs_status = EXIT_FAILURE;
     }
diff --git a/Applications/CLI/ogs_python_bindings.cpp b/Applications/CLI/ogs_python_bindings.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..12c0225bf053f874acdb93390faa48bf6da2482e
--- /dev/null
+++ b/Applications/CLI/ogs_python_bindings.cpp
@@ -0,0 +1,38 @@
+#ifdef OGS_USE_PYTHON
+
+#include <pybind11/embed.h>
+
+#include "ogs_python_bindings.h"
+
+#ifndef OGS_BUILD_SHARED_LIBS
+extern "C" PyObject* pybind11_init_impl_OpenGeoSys();
+
+// Hackish trick that hopefully ensures that the linker won't strip the symbol
+// pointed to by p from the library being built.
+template <typename T>
+void mark_used(T p)
+{
+    volatile T vp = (volatile T)p;
+    vp = vp;
+}
+
+#endif  // OGS_BUILD_SHARED_LIBS
+
+namespace ApplicationsLib
+{
+pybind11::scoped_interpreter setupEmbeddedPython()
+{
+#ifndef OGS_BUILD_SHARED_LIBS
+    // pybind11_init_impl_OpenGeoSys is the function initializing the embedded
+    // OpenGeoSys Python module. The name is generated by pybind11. If it is not
+    // obvious that this symbol is actually used, the linker might remove it
+    // under certain circumstances.
+    mark_used(&pybind11_init_impl_OpenGeoSys);
+#endif
+
+    return pybind11::scoped_interpreter{};
+}
+
+}  // namespace ApplicationsLib
+
+#endif  // OGS_USE_PYTHON
diff --git a/Applications/CLI/ogs_python_bindings.h b/Applications/CLI/ogs_python_bindings.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4b1e2a9442bc55800a9a5b5d8410f1ab2428279
--- /dev/null
+++ b/Applications/CLI/ogs_python_bindings.h
@@ -0,0 +1,15 @@
+#pragma once
+
+#ifdef OGS_USE_PYTHON
+
+#include <pybind11/embed.h>
+
+namespace ApplicationsLib
+{
+/// Sets up an embedded Python interpreter and makes sure that the OpenGeoSys
+/// Python module is not removed by the linker.
+pybind11::scoped_interpreter setupEmbeddedPython();
+
+}  // namespace ApplicationsLib
+
+#endif  // OGS_USE_PYTHON
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3d87e8ba12738e631d4dfbb34c53744ff2ab4dc4..d183b58ad142567857f6cc3daa539fcaa6f691f6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -210,6 +210,8 @@ option(OGS_ENABLE_ELEMENT_PYRAMID "Build FEM elements for pyramids." ON)
 
 option(OGS_CHECK_HEADER_COMPILATION "Check header for standalone compilation." OFF)
 
+option(OGS_USE_PYTHON "Interface with python" OFF)
+
 ###################
 ### Definitions ###
 ###################
@@ -259,12 +261,16 @@ if(OGS_USE_EIGEN)
         add_definitions(-DEIGEN_INITIALIZE_MATRICES_BY_NAN)
     endif()
 endif()
-
 if(MSVC AND OGS_32_BIT)
     add_definitions(-DEIGEN_MAX_ALIGN_BYTES=0 -DEIGEN_DONT_ALIGN)
 endif()
 # End Eigen
 
+if(OGS_USE_PYTHON)
+    add_definitions(-DOGS_USE_PYTHON)
+    find_package(pybind11 REQUIRED)
+endif()
+
 if (OGS_FATAL_ABORT)
     add_definitions(-DOGS_FATAL_ABORT)
 endif()