diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_geometry.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_geometry.md
index 7bb53b68eb98d1a010e70be0a362ed1ef025fc76..5ded223a30d2a226a5060bcead103d46059e356d 100644
--- a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_geometry.md
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/t_geometry.md
@@ -1 +1 @@
-The name of the on which the boundary condition is defined.
+The name of the geometry on which the boundary condition is defined.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/i_source_terms.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/i_source_terms.md
new file mode 100644
index 0000000000000000000000000000000000000000..8c5caaa4a7e38517d9438f237505661ace9ae7e3
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/i_source_terms.md
@@ -0,0 +1 @@
+Contains source terms for the process variable.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/c_Nodal.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/c_Nodal.md
new file mode 100644
index 0000000000000000000000000000000000000000..a4f9ba28692f8a0152fdebbef12dcd6e64550aa6
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/c_Nodal.md
@@ -0,0 +1 @@
+Declares a source term to be associated to a node.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/t_value.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/t_value.md
new file mode 100644
index 0000000000000000000000000000000000000000..c27bd71e47f57aed6873e4416247e165e63564d7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Nodal/t_value.md
@@ -0,0 +1 @@
+The value for the nodal source term.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/i_source_term.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/i_source_term.md
new file mode 100644
index 0000000000000000000000000000000000000000..60619b9e1c089c5aa742046024ddd226444160e1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/i_source_term.md
@@ -0,0 +1 @@
+Defines a source term for a process variable.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_component.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_component.md
new file mode 100644
index 0000000000000000000000000000000000000000..71df2b71742809e65402d4b3ef37a0fc7050cba3
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_component.md
@@ -0,0 +1,3 @@
+For the multi-component process variables specifies the component of the
+source term. The component ids start with 0. For single-component
+variables this is optional and defaults to 0.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometrical_set.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometrical_set.md
new file mode 100644
index 0000000000000000000000000000000000000000..66588b680141f84325ac570607e76fff6b1c9087
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometrical_set.md
@@ -0,0 +1,4 @@
+The name of the geometrical set (see \ref ogs_file_param__gml__name) in which
+the \ref
+ogs_file_param__prj__process_variables__process_variable__source_terms__source_term__geometry
+"geometry" is defined.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometry.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometry.md
new file mode 100644
index 0000000000000000000000000000000000000000..b792a76905a9661d786d7c03b7708bff79d722e4
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_geometry.md
@@ -0,0 +1 @@
+The name of the geometry on which the source term is defined.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_type.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_type.md
new file mode 100644
index 0000000000000000000000000000000000000000..71e7165a64cbbb6610b3c8a2f42b23dc495f30ab
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/t_type.md
@@ -0,0 +1 @@
+Sets the type of the source term.
diff --git a/MeshGeoToolsLib/CreateSearchLength.cpp b/MeshGeoToolsLib/CreateSearchLength.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b10264e8d0fd4d034d68020bfb9d25ef6575d6bc
--- /dev/null
+++ b/MeshGeoToolsLib/CreateSearchLength.cpp
@@ -0,0 +1,49 @@
+/**
+ * @file
+ * @copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "CreateSearchLength.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "MeshGeoToolsLib/HeuristicSearchLength.h"
+#include "MeshGeoToolsLib/SearchLength.h"
+
+namespace MeshGeoToolsLib
+{
+std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
+    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh)
+{
+    boost::optional<BaseLib::ConfigTree> config =
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm}
+        external_config.getConfigSubtreeOptional("search_length_algorithm");
+
+    if (!config)
+        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
+            new MeshGeoToolsLib::SearchLength()};
+
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__type}
+    std::string const type = config->getConfigParameter<std::string>("type");
+
+    if (type == "fixed")
+    {
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__fixed__value}
+        double const length = config->getConfigParameter<double>("value");
+        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
+            new MeshGeoToolsLib::SearchLength(length)};
+    }
+    if (type == "heuristic")
+    {
+        return std::unique_ptr<MeshGeoToolsLib::HeuristicSearchLength>{
+            new MeshGeoToolsLib::HeuristicSearchLength(mesh)};
+    }
+    OGS_FATAL("Unknown search length algorithm type '%s'.", type.c_str());
+}
+
+}  // end namespace MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/CreateSearchLength.h b/MeshGeoToolsLib/CreateSearchLength.h
index a46a01baa90ba7af411b2fff4d84c017c6b76a5b..475f71b94be33db565e2c94353d50e0d4cfc8a4b 100644
--- a/MeshGeoToolsLib/CreateSearchLength.h
+++ b/MeshGeoToolsLib/CreateSearchLength.h
@@ -13,40 +13,20 @@
 
 #include <memory>
 
-#include "BaseLib/ConfigTree.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
 
-#include "MeshGeoToolsLib/HeuristicSearchLength.h"
-#include "MeshGeoToolsLib/SearchLength.h"
+namespace MeshLib
+{
+class Mesh;
+}
 
 namespace MeshGeoToolsLib
 {
-std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
-    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh)
-{
-    boost::optional<BaseLib::ConfigTree> config =
-        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm}
-        external_config.getConfigSubtreeOptional("search_length_algorithm");
-
-    if (!config)
-        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
-            new MeshGeoToolsLib::SearchLength()};
-
-    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__type}
-    std::string const type = config->getConfigParameter<std::string>("type");
-
-    if (type == "fixed")
-    {
-        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__search_length_algorithm__fixed__value}
-        double const length = config->getConfigParameter<double>("value");
-        return std::unique_ptr<MeshGeoToolsLib::SearchLength>{
-            new MeshGeoToolsLib::SearchLength(length)};
-    }
-    if (type == "heuristic")
-    {
-        return std::unique_ptr<MeshGeoToolsLib::HeuristicSearchLength>{
-            new MeshGeoToolsLib::HeuristicSearchLength(mesh)};
-    }
-    OGS_FATAL("Unknown search length algorithm type '%s'.", type.c_str());
-}
+class SearchLength;
 
+std::unique_ptr<MeshGeoToolsLib::SearchLength> createSearchLengthAlgorithm(
+    BaseLib::ConfigTree const& external_config, MeshLib::Mesh const& mesh);
 }  // end namespace MeshGeoToolsLib
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 75ec2adcc563fee3e91ac42a3a5796e6787a1216..6fd7c6b13dc9cc18d62978452541a848360dfc9a 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -21,6 +21,7 @@ APPEND_SOURCE_FILES(SOURCES Parameter)
 APPEND_SOURCE_FILES(SOURCES PhaseField)
 APPEND_SOURCE_FILES(SOURCES RichardsFlow)
 APPEND_SOURCE_FILES(SOURCES SmallDeformation)
+APPEND_SOURCE_FILES(SOURCES SourceTerms)
 APPEND_SOURCE_FILES(SOURCES TES)
 APPEND_SOURCE_FILES(SOURCES ThermoMechanics)
 APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPP)
diff --git a/ProcessLib/GroundwaterFlow/Tests.cmake b/ProcessLib/GroundwaterFlow/Tests.cmake
index ff7582f71d535164c9fc15520b423b71f819a433..3b06c66ad4c9c55bf9a0383d2f01ae86fb1219e4 100644
--- a/ProcessLib/GroundwaterFlow/Tests.cmake
+++ b/ProcessLib/GroundwaterFlow/Tests.cmake
@@ -597,3 +597,106 @@ AddTest(
     DIFF_DATA
     expected_dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu pressure pressure
 )
+
+# tests for nodal source term implementation
+# For the special setup with a dirac source term at position (xi, eta) the
+# analytical solution in 2 dimensions is valid:
+# u(x,y) = ln(sqrt((x-xi)^2+(y-eta)^2))/(2 * Pi)
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e1
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e1_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 0.7 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e1_expected.vtu circle_1e1_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e1_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e2
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e2_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.1 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e2_expected.vtu circle_1e2_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e2_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e3
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e3_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.6 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e3_expected.vtu circle_1e3_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e3_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NodalSourceTerm_circle_1e4
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e4_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.8 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e4_expected.vtu circle_1e4_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e4_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_circle_1e5
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e5_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 2.15 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e5_expected.vtu circle_1e5_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e5_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_circle_1e6
+    PATH Elliptic/circle_radius_1
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS circle_1e6_axi.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 2.52 RELTOL 1e-16
+    DIFF_DATA
+    line_1_lines_1e6_expected.vtu circle_1e6_axi_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS circle_1e6_axi_pcs_0_ts_1_t_1.000000.vtu
+)
+
+AddTest(
+    NAME LARGE_GroundWaterFlowProcess_NodalSourceTerm_square_1e6
+    PATH Elliptic/square_1x1_GroundWaterFlow
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e6_with_nodal_sources.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    ABSTOL 1.41 RELTOL 1e-16
+    DIFF_DATA
+    square_1x1_quad_1e6_nodal_sources_expected.vtu square_1e6_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure
+    VIS square_1e6_pcs_0_ts_1_t_1.000000.vtu
+)
+
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index afe758b7ebb500884aa67674dcffebbbb75dd260..9d066191a12a2b8d1d095835d88f3424068fad0d 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -59,6 +59,19 @@ void Process::initialize()
     DBUG("Initialize boundary conditions.");
     _boundary_conditions.addBCsForProcessVariables(
         _process_variables, *_local_to_global_index_map, _integration_order);
+
+    for (int variable_id = 0;
+         variable_id < static_cast<int>(_process_variables.size());
+         ++variable_id)
+    {
+        ProcessVariable& pv = _process_variables[variable_id];
+        auto sts =
+            pv.createSourceTerms(*_local_to_global_index_map, variable_id,
+                                 _integration_order);
+
+        std::move(sts.begin(), sts.end(),
+                  std::back_inserter(_source_terms));
+    }
 }
 
 void Process::setInitialConditions(double const t, GlobalVector& x)
@@ -128,6 +141,11 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
     assembleConcreteProcess(t, x, M, K, b);
 
     _boundary_conditions.applyNaturalBC(t, x, K, b);
+
+    for (auto const& st : _source_terms)
+    {
+        st->integrateNodalSourceTerm(t, b);
+    }
 }
 
 void Process::assembleWithJacobian(const double t, GlobalVector const& x,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 88ef33bf5ee60e203fba6c30ad7822211c3ec475..b9fba6a81f6c91bb7e5d1fcbafa1d0274f0fb353 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -14,6 +14,7 @@
 #include "NumLib/ODESolver/TimeDiscretization.h"
 #include "NumLib/NamedFunctionCaller.h"
 #include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h"
+#include "ProcessLib/SourceTerms/NodalSourceTerm.h"
 #include "ProcessLib/Parameter/Parameter.h"
 
 #include "ExtrapolatorData.h"
@@ -218,6 +219,7 @@ private:
     std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
 
     BoundaryConditionCollection _boundary_conditions;
+    std::vector<std::unique_ptr<NodalSourceTerm>> _source_terms;
 
     ExtrapolatorData _extrapolator_data;
 };
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 404c978e2dfda928bf9e87558f743c1c89e6b879..987e6b0989f3659a7254b87ac36044769de77230 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -33,7 +33,8 @@ ProcessVariable::ProcessVariable(
           //! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
           config.getConfigParameter<std::string>("initial_condition"),
           parameters, _n_components)),
-      _bc_builder(std::make_unique<BoundaryConditionBuilder>())
+      _bc_builder(std::make_unique<BoundaryConditionBuilder>()),
+      _source_term_builder(std::make_unique<SourceTermBuilder>())
 {
     DBUG("Constructing process variable %s", _name.c_str());
 
@@ -82,6 +83,49 @@ ProcessVariable::ProcessVariable(
     } else {
         INFO("No boundary conditions found.");
     }
+
+    // Source terms
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms}
+    if (auto sts_config = config.getConfigSubtreeOptional("source_terms"))
+    {
+        for (auto st_config :
+             //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term}
+             sts_config->getConfigSubtreeList("source_term"))
+        {
+            auto const geometrical_set_name =
+                    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometrical_set}
+                   st_config.getConfigParameter<std::string>("geometrical_set");
+            auto const geometry_name =
+                    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometry}
+                    st_config.getConfigParameter<std::string>("geometry");
+
+            GeoLib::GeoObject const* const geometry =
+                geometries.getGeoObject(geometrical_set_name, geometry_name);
+
+            if (! geometry)
+                OGS_FATAL(
+                    "No geometry with name `%s' has been found in the "
+                    "geometrical set `%s'.",
+                    geometry_name.c_str(), geometrical_set_name.c_str());
+
+            DBUG(
+                "Found geometry type \"%s\"",
+                GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str());
+
+            auto component_id =
+                //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__component}
+                st_config.getConfigParameterOptional<int>("component");
+
+            if (!component_id && _n_components == 1)
+                // default value for single component vars.
+                component_id = 0;
+
+            _source_term_configs.emplace_back(std::move(st_config), *geometry,
+                                              component_id);
+        }
+    } else {
+        INFO("No source terms found.");
+    }
 }
 
 ProcessVariable::ProcessVariable(ProcessVariable&& other)
@@ -91,7 +135,9 @@ ProcessVariable::ProcessVariable(ProcessVariable&& other)
       _shapefunction_order(other._shapefunction_order),
       _initial_condition(std::move(other._initial_condition)),
       _bc_configs(std::move(other._bc_configs)),
-      _bc_builder(std::move(other._bc_builder))
+      _bc_builder(std::move(other._bc_builder)),
+      _source_term_configs(std::move(other._source_term_configs)),
+      _source_term_builder(std::move(other._source_term_builder))
 {
 }
 
@@ -128,4 +174,20 @@ ProcessVariable::createBoundaryConditions(
     return bcs;
 }
 
+std::vector<std::unique_ptr<NodalSourceTerm>>
+ProcessVariable::createSourceTerms(
+    const NumLib::LocalToGlobalIndexMap& dof_table,
+    const int variable_id,
+    unsigned const integration_order)
+{
+    std::vector<std::unique_ptr<NodalSourceTerm>> source_terms;
+
+    for (auto& config : _source_term_configs)
+        source_terms.emplace_back(_source_term_builder->createSourceTerm(
+            config, dof_table, _mesh, variable_id, integration_order,
+            _shapefunction_order));
+
+    return source_terms;
+}
+
 }  // namespace ProcessLib
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 56aa3c9e797d31bd4e4223d4c1c417e4d2d45194..1fe809b6c45a9eaabdc560d8fe85721cd7487ee8 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -12,6 +12,8 @@
 #include "ProcessLib/BoundaryCondition/BoundaryCondition.h"
 #include "ProcessLib/BoundaryCondition/BoundaryConditionConfig.h"
 #include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/SourceTerms/SourceTermConfig.h"
+#include "ProcessLib/SourceTerms/SourceTermBuilder.h"
 
 namespace MeshLib
 {
@@ -53,6 +55,10 @@ public:
         unsigned const integration_order,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters);
 
+    std::vector<std::unique_ptr<NodalSourceTerm>> createSourceTerms(
+        const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id,
+        unsigned const integration_order);
+
     Parameter<double> const& getInitialCondition() const
     {
         return _initial_condition;
@@ -87,6 +93,8 @@ private:
 
     std::vector<BoundaryConditionConfig> _bc_configs;
     std::unique_ptr<BoundaryConditionBuilder> _bc_builder;
+    std::vector<SourceTermConfig> _source_term_configs;
+    std::unique_ptr<SourceTermBuilder> _source_term_builder;
 };
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp b/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..990703d901467e8b46e9378cf777ebe3a70bab64
--- /dev/null
+++ b/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp
@@ -0,0 +1,35 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateNodalSourceTerm.h"
+
+#include "BaseLib/FileTools.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "NodalSourceTerm.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NodalSourceTerm> createNodalSourceTerm(
+    BaseLib::ConfigTree const& config,
+    const NumLib::LocalToGlobalIndexMap& dof_table, std::size_t const mesh_id,
+    std::size_t const node_id, const int variable_id, const int component_id)
+{
+    DBUG("Constructing NodalSourceTerm from config.");
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+    config.checkConfigParameter("type", "Nodal");
+
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__Nodal__value}
+    auto const nodal_value = config.getConfigParameter<double>("value");
+    DBUG("Using value %f as nodal source term", nodal_value);
+
+    return std::make_unique<NodalSourceTerm>(
+        dof_table, mesh_id, node_id, variable_id, component_id, nodal_value);
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/CreateNodalSourceTerm.h b/ProcessLib/SourceTerms/CreateNodalSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..e480133b5d6a508d28052a844695c8295dd7112f
--- /dev/null
+++ b/ProcessLib/SourceTerms/CreateNodalSourceTerm.h
@@ -0,0 +1,22 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NodalSourceTerm> createNodalSourceTerm(
+    BaseLib::ConfigTree const& config,
+    const NumLib::LocalToGlobalIndexMap& dof_table, std::size_t mesh_id,
+    std::size_t const node_id, const int variable_id, const int component_id);
+
+}   // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.cpp b/ProcessLib/SourceTerms/NodalSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..874d6a136f6e8730be82055575a07686ef3f97fd
--- /dev/null
+++ b/ProcessLib/SourceTerms/NodalSourceTerm.cpp
@@ -0,0 +1,43 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "NodalSourceTerm.h"
+
+#include <cassert>
+
+namespace ProcessLib
+{
+NodalSourceTerm::NodalSourceTerm(const NumLib::LocalToGlobalIndexMap& dof_table,
+                                 std::size_t const mesh_id,
+                                 std::size_t const node_id,
+                                 const int variable_id, const int component_id,
+                                 double value)
+    : _dof_table(dof_table),
+      _mesh_id(mesh_id),
+      _node_id(node_id),
+      _variable_id(variable_id),
+      _component_id(component_id),
+      _value(value)
+{
+    DBUG("Create NodalSourceTerm.");
+}
+
+void NodalSourceTerm::integrateNodalSourceTerm(
+    const double t,
+    GlobalVector& b) const
+{
+    DBUG("Assemble NodalSourceTerm.");
+
+    MeshLib::Location const l{_mesh_id, MeshLib::MeshItemType::Node, _node_id};
+    auto const index =
+        _dof_table.getGlobalIndex(l, _variable_id, _component_id);
+    b.add(index, _value);
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.h b/ProcessLib/SourceTerms/NodalSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4a7ef4622fb234ae874bd727c28e33d7a1b905f
--- /dev/null
+++ b/ProcessLib/SourceTerms/NodalSourceTerm.h
@@ -0,0 +1,37 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace ProcessLib
+{
+class NodalSourceTerm final
+{
+public:
+    NodalSourceTerm(const NumLib::LocalToGlobalIndexMap& dof_table,
+                    std::size_t const mesh_id, std::size_t const node_id,
+                    const int variable_id, const int component_id,
+                    double value);
+
+    void integrateNodalSourceTerm(
+        const double t,
+        GlobalVector& b) const;
+
+private:
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+    std::size_t const _mesh_id;
+    std::size_t const _node_id;
+    int const _variable_id;
+    int const _component_id;
+    double const _value;
+};
+
+}   // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/SourceTermBuilder.cpp b/ProcessLib/SourceTerms/SourceTermBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a9340858c0d2843e4fe0885790eaf4d92f74917
--- /dev/null
+++ b/ProcessLib/SourceTerms/SourceTermBuilder.cpp
@@ -0,0 +1,98 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "SourceTermBuilder.h"
+#include "SourceTermConfig.h"
+#include "CreateNodalSourceTerm.h"
+#include "NodalSourceTerm.h"
+#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshGeoToolsLib/CreateSearchLength.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshGeoToolsLib/SearchLength.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NodalSourceTerm> SourceTermBuilder::createSourceTerm(
+    const SourceTermConfig& config,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+    const int variable_id, const unsigned integration_order,
+    const unsigned shapefunction_order)
+{
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+    auto const type = config.config.peekConfigParameter<std::string>("type");
+
+    if (type == "Nodal")
+    {
+        return createNodalSourceTerm(config, dof_table, mesh, variable_id,
+                                     integration_order, shapefunction_order);
+    }
+
+    OGS_FATAL("Unknown source term type: `%s'.", type.c_str());
+}
+
+std::unique_ptr<NodalSourceTerm> SourceTermBuilder::createNodalSourceTerm(
+    const SourceTermConfig& config,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+    const int variable_id, const unsigned /*integration_order*/,
+    const unsigned /*shapefunction_order*/)
+{
+    std::unique_ptr<MeshGeoToolsLib::SearchLength> search_length_algorithm =
+        MeshGeoToolsLib::createSearchLengthAlgorithm(config.config, mesh);
+
+    MeshGeoToolsLib::MeshNodeSearcher const& mesh_node_searcher =
+        MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+            mesh, std::move(search_length_algorithm));
+
+    // Find nodes' ids on the given mesh on which this source term is defined.
+    std::vector<std::size_t> ids =
+        mesh_node_searcher.getMeshNodeIDs(config.geometry);
+
+    // Filter out ids, which are not part of mesh subsets corresponding to
+    // the variable_id and component_id.
+
+    // Sorted ids of all mesh_subsets.
+    std::vector<std::size_t> sorted_nodes_ids;
+
+    auto const& mesh_subsets =
+        dof_table.getMeshSubsets(variable_id, *config.component_id);
+    for (auto const& mesh_subset : mesh_subsets)
+    {
+        auto const& nodes = mesh_subset->getNodes();
+        sorted_nodes_ids.reserve(sorted_nodes_ids.size() + nodes.size());
+        std::transform(std::begin(nodes), std::end(nodes),
+                       std::back_inserter(sorted_nodes_ids),
+                       [](MeshLib::Node* const n) { return n->getID(); });
+    }
+    std::sort(std::begin(sorted_nodes_ids), std::end(sorted_nodes_ids));
+
+    auto ids_new_end_iterator = std::end(ids);
+    ids_new_end_iterator = std::remove_if(
+        std::begin(ids), ids_new_end_iterator,
+        [&sorted_nodes_ids](std::size_t const node_id) {
+            return !std::binary_search(std::begin(sorted_nodes_ids),
+                                       std::end(sorted_nodes_ids), node_id);
+        });
+    ids.erase(ids_new_end_iterator, std::end(ids));
+
+    DBUG(
+        "Found %d nodes for nodal source term for the variable %d and "
+        "component %d",
+        ids.size(), variable_id, *config.component_id);
+
+    if (ids.size() != 1)
+        OGS_FATAL(
+            "Found %d nodes for nodal source term, but exactly one node is "
+            "required.");
+
+    return ProcessLib::createNodalSourceTerm(
+        config.config, dof_table, mesh.getID(), ids[0], variable_id,
+        *config.component_id);
+}
+
+}  // ProcessLib
diff --git a/ProcessLib/SourceTerms/SourceTermBuilder.h b/ProcessLib/SourceTerms/SourceTermBuilder.h
new file mode 100644
index 0000000000000000000000000000000000000000..dbc0981f8f160af066cc662ec92fd6d92b324804
--- /dev/null
+++ b/ProcessLib/SourceTerms/SourceTermBuilder.h
@@ -0,0 +1,60 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NodalSourceTerm.h"
+
+namespace GeoLib
+{
+class GeoObject;
+}
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+}
+
+namespace MeshGeoToolsLib
+{
+class BoundaryElementsSearcher;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+template <typename>
+struct IndexValueVector;
+}
+
+namespace ProcessLib
+{
+struct SourceTermConfig;
+
+class SourceTermBuilder
+{
+public:
+    virtual ~SourceTermBuilder() = default;
+
+    virtual std::unique_ptr<NodalSourceTerm> createSourceTerm(
+        const SourceTermConfig& config,
+        const NumLib::LocalToGlobalIndexMap& dof_table,
+        const MeshLib::Mesh& mesh, const int variable_id,
+        const unsigned integration_order, const unsigned shapefunction_order);
+
+protected:
+    virtual std::unique_ptr<NodalSourceTerm> createNodalSourceTerm(
+        const SourceTermConfig& config,
+        const NumLib::LocalToGlobalIndexMap& dof_table,
+        const MeshLib::Mesh& mesh, const int variable_id,
+        const unsigned integration_order, const unsigned shapefunction_order);
+};
+
+}  // ProcessLib
diff --git a/ProcessLib/SourceTerms/SourceTermConfig.h b/ProcessLib/SourceTerms/SourceTermConfig.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d495a6d2b9048cf5308563947c528aa4483ff14
--- /dev/null
+++ b/ProcessLib/SourceTerms/SourceTermConfig.h
@@ -0,0 +1,40 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "GeoLib/GEOObjects.h"
+
+namespace ProcessLib
+{
+struct SourceTermConfig final
+{
+    SourceTermConfig(BaseLib::ConfigTree&& config_,
+                     GeoLib::GeoObject const& geometry_,
+                     boost::optional<int> const component_id_)
+        : config(std::move(config_)),
+          geometry(geometry_),
+          component_id(component_id_)
+    {
+    }
+
+    SourceTermConfig(SourceTermConfig&& other)
+        : config(std::move(other.config)),
+          geometry(other.geometry),
+          component_id(other.component_id)
+    {
+    }
+
+    BaseLib::ConfigTree config;
+    GeoLib::GeoObject const& geometry;
+    boost::optional<int> const component_id;
+};
+
+}  // ProcessLib
diff --git a/Tests/Data b/Tests/Data
index a9f0365e2ac73714892c4b6c103af93b4141d8bf..8a91a43ac885fbcb3876a808708ec588b9fea42c 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit a9f0365e2ac73714892c4b6c103af93b4141d8bf
+Subproject commit 8a91a43ac885fbcb3876a808708ec588b9fea42c
diff --git a/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_analytical_solution_head.png b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_analytical_solution_head.png
new file mode 100644
index 0000000000000000000000000000000000000000..77f3e513bcc2433250fd8888d7437a7c05a40278
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_analytical_solution_head.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9081be9dbb3bc00928af83861c0d5899db6ae2568986eb4c4e5b7b7117459523
+size 12958
diff --git a/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head.png b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head.png
new file mode 100644
index 0000000000000000000000000000000000000000..ea2902842748f5c7b066518238c513d407043993
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5a50ac5f92ebad2ceb3c1f61ceb9158e5f90c89e79eab4f1072be00d61d5a3ac
+size 13548
diff --git a/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head_log_scale.png b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head_log_scale.png
new file mode 100644
index 0000000000000000000000000000000000000000..95fe75a1a667f04391c73de2212091fa8e31d1aa
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head_log_scale.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bdff20c75f67452f136d01f56ff25bc1c96d8d4a5c2aa783e720421be4170159
+size 20683
diff --git a/web/content/docs/benchmarks/elliptic/groundwater-flow-dirichlet-nodal-source-term.md b/web/content/docs/benchmarks/elliptic/groundwater-flow-dirichlet-nodal-source-term.md
new file mode 100644
index 0000000000000000000000000000000000000000..1e06da84eaae98746a299015142fc4d736da187b
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/groundwater-flow-dirichlet-nodal-source-term.md
@@ -0,0 +1,82 @@
++++
+date = "2017-02-15T11:17:39+01:00"
+title = "Groundwater Flow (Nodal Source Term)"
+project = "Elliptic/circle_radius_1/circle_1e6_axi.prj"
+author = "Thomas Fischer"
+weight = 101
+
+aliases = [ "/docs/benchmarks/" ] # First benchmark page
+
+[menu]
+  [menu.benchmarks]
+    parent = "elliptic"
+
++++
+
+{{< project-link >}}
+
+## Equations
+
+We solve the Poisson equation:
+$$
+\begin{equation}
+k\; \Delta h = f(x) \quad \text{in }\Omega
+\end{equation}$$
+w.r.t boundary conditions
+$$
+\eqalign{
+h(x) = g_D(x) &\quad \text{on }\Gamma_D,\cr
+k\;{\partial h(x) \over \partial n} = g_N(x) &\quad \text{on }\Gamma_N,
+}$$
+
+where $h$ could be hydraulic head, the subscripts $D$ and $N$ denote the Dirichlet- and Neumann-type boundary conditions, $n$ is the normal vector pointing outside of $\Omega$, and $\Gamma = \Gamma_D \cup \Gamma_N$ and $\Gamma_D \cap \Gamma_N = \emptyset$.
+
+## Problem specification and analytical solution
+
+We solve the Poisson equation on a circle domain with radius $r = 1$ with $k = 1$ w.r.t. the specific boundary conditions:
+$$
+\eqalign{
+h(x,y) = 0 &\quad \text{on } (x^2 + y^2 = 1) \subset \Gamma_D,\cr
+}$$
+The solution of this problem is
+$$
+h(x,y) = \int \int f(\xi, \eta) G(x, y) d \xi d \eta,
+$$
+where $G(x, y)$ is the Green's function. For the example at hand $G(x, y)$ is:
+$$
+G(x, y) = \frac{1}{2 \pi} \ln \sqrt{(x-\xi)^2 + (y-\eta)^2}.
+$$
+With a nodal source term of 1 at $(0.0, 0.0)$ the analytical solution is
+$$
+h(x,y) = \frac{1}{2 \pi} \ln \sqrt{x^2 + y^2}.
+$$
+
+
+## Input files
+
+The main project file is `square_1e6_with_nodal_sources.prj`. It describes the process to be solved and the related process variables together with their initial and boundary conditions as well as the definition of the nodal source term. It also references the mesh and geometrical objects defined on the mesh.
+
+The geometries used to specify the boundary conditions and the source term are given in the `square_1x1.gml` file.
+
+The input mesh `square_1x1_quad_1e6.vtu` is stored in the VTK file format and can be directly visualized in Paraview for example.
+
+## Running simulation
+
+To start the simulation (after successful compilation) run:
+```bash
+$ ogs circle_1e6_axi.prj
+```
+
+It will produce some output and write the computed result into a data array of the written vtu file.
+
+
+## Results and evaluation
+
+### Comparison of the analytical solution and the computed solution
+
+{{< img src="../circle_1e6_gwf_with_nodal_source_term_analytical_solution_head.png" >}}
+
+{{< img src="../circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head.png" >}}
+
+{{< img src="../circle_1e6_gwf_with_nodal_source_term_diff_analytical_solution_head_log_scale.png" >}}
+