diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index 11034eb047086d03c50107e8f7756fe83d0646bd..08be63a171381d0cd44382213aa281e8bd50d6fc 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -150,33 +150,6 @@ public:
                 this->_integration_order);
     }
 
-    void initializeBoundaryConditions() override
-    {
-        MeshGeoToolsLib::MeshNodeSearcher& hydraulic_head_mesh_node_searcher =
-            MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
-                _hydraulic_head->getMesh());
-
-        //
-        // Neumann boundary conditions.
-        //
-        {
-            // Find mesh nodes.
-            MeshGeoToolsLib::BoundaryElementsSearcher hydraulic_head_mesh_element_searcher(
-                _hydraulic_head->getMesh(), hydraulic_head_mesh_node_searcher);
-
-            // Create a neumann BC for the hydraulic head storing them in the
-            // _neumann_bcs vector.
-            _hydraulic_head->createNeumannBcs(
-                    std::back_inserter(this->_neumann_bcs),
-                    hydraulic_head_mesh_element_searcher,
-                    this->_global_setup,
-                    this->_integration_order,
-                    *this->_local_to_global_index_map,
-                    0,
-                    *this->_mesh_subset_all_nodes);
-        }
-    }
-
     void initializeMeshSubsets(MeshLib::Mesh const& mesh) override
     {
         // Create single component dof in every of the mesh's nodes.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index ec30aea641e3e8ca321b95a0ab742f9816c1c5e8..cd4f87884dd2f8a687ec02a0050077a94c2c25b5 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -63,8 +63,6 @@ public:
 	/// Creates mesh subsets, i.e. components, for given mesh.
 	virtual void initializeMeshSubsets(MeshLib::Mesh const& mesh) = 0;
 
-	virtual void initializeBoundaryConditions() = 0;
-
 	void initialize()
 	{
 		DBUG("Initialize process.");
@@ -95,10 +93,9 @@ public:
 			DBUG("Initialize boundary conditions.");
 			createDirichletBcs(*pv, 0);  // 0 is the component id
 
+			createNeumannBcs(*pv, 0);  // 0 is the component id
 		}
 
-		initializeBoundaryConditions();
-
 		for (auto& bc : _neumann_bcs)
 			bc->initialize(_global_setup, *_A, *_rhs, _mesh.getDimension());
 	}
@@ -161,6 +158,26 @@ protected:
 		                                 component_id);
 	}
 
+	void createNeumannBcs(ProcessVariable& variable, int const component_id)
+	{
+		// Find mesh nodes.
+		MeshGeoToolsLib::MeshNodeSearcher& mesh_node_searcher =
+		    MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+		        variable.getMesh());
+		MeshGeoToolsLib::BoundaryElementsSearcher mesh_element_searcher(
+		    variable.getMesh(), mesh_node_searcher);
+
+		// Create a neumann BC for the process variable storing them in the
+		// _neumann_bcs vector.
+		variable.createNeumannBcs(std::back_inserter(_neumann_bcs),
+		                          mesh_element_searcher,
+		                          _global_setup,
+		                          _integration_order,
+		                          *_local_to_global_index_map,
+		                          component_id,
+		                          *_mesh_subset_all_nodes);
+	}
+
 private:
 	/// Creates global matrix, rhs and solution vectors, and the linear solver.
 	void createLinearSolver(std::string const& solver_name)