diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index 462dc8acd37a5159c62f4b3ddf9ad573e84a2a33..04d30ad53718fdae47409b46ec5dc8a059a30ddd 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -39,6 +39,7 @@
 #include "BaseLib/ConfigTreeUtil.h"
 #include "BaseLib/DateTools.h"
 #include "BaseLib/FileTools.h"
+#include "BaseLib/PrjProcessing.h"
 #include "BaseLib/RunTime.h"
 #include "InfoLib/CMakeInfo.h"
 #include "InfoLib/GitInfo.h"
@@ -89,6 +90,12 @@ int main(int argc, char* argv[])
         false, "");
     cmd.add(xml_patch_files);
 
+    TCLAP::SwitchArg write_prj("",
+                               "write-prj",
+                               "Writes processed project file to output "
+                               "path / [prj_base_name]_processed.prj.");
+    cmd.add(write_prj);
+
     TCLAP::ValueArg<std::string> outdir_arg("o", "output-directory",
                                             "the output directory to write to",
                                             false, "", "PATH");
@@ -195,9 +202,14 @@ int main(int argc, char* argv[])
 
             run_time.start();
 
+            std::stringstream prj_stream;
+            BaseLib::prepareProjectFile(
+                prj_stream, project_arg.getValue(), xml_patch_files.getValue(),
+                write_prj.getValue(), outdir_arg.getValue());
+
             auto project_config = BaseLib::makeConfigTree(
                 project_arg.getValue(), !nonfatal_arg.getValue(),
-                "OpenGeoSysProject", xml_patch_files.getValue());
+                "OpenGeoSysProject", prj_stream);
 
             BaseLib::setProjectDirectory(
                 BaseLib::extractPath(project_arg.getValue()));
diff --git a/BaseLib/ConfigTreeUtil.cpp b/BaseLib/ConfigTreeUtil.cpp
index 954fc5189cb4d2d4332325581524b1ea299236bb..786692aff08d9ce09e33861f796c6f9cef9eda36 100644
--- a/BaseLib/ConfigTreeUtil.cpp
+++ b/BaseLib/ConfigTreeUtil.cpp
@@ -10,220 +10,18 @@
 
 #include "ConfigTreeUtil.h"
 
-#include <libxml/parser.h>
-#include <libxml/tree.h>
-#include <spdlog/spdlog.h>
-#include <xml_patch.h>
-
 #include <boost/property_tree/xml_parser.hpp>
-#include <filesystem>
 #include <regex>
 
-#include "BaseLib/FileTools.h"
 #include "Error.h"
 #include "Logging.h"
 
 namespace BaseLib
 {
-// Adapted from
-// https://stackoverflow.com/questions/8154107/how-do-i-merge-update-a-boostproperty-treeptree/8175833
-template <typename T>
-void traverse_recursive(
-    boost::property_tree::ptree& parent,
-    boost::property_tree::ptree::path_type const& child_path,
-    boost::property_tree::ptree& child,
-    std::filesystem::path const& bench_dir,
-    T& method)
-{
-    using boost::property_tree::ptree;
-
-    method(parent, child_path, child, bench_dir);
-    for (auto& [key, tree] : child)
-    {
-        ptree::path_type const cur_path = child_path / ptree::path_type(key);
-        traverse_recursive(child, cur_path, tree, bench_dir, method);
-    }
-}
-
-template <typename T>
-void traverse(boost::property_tree::ptree& parent,
-              const std::filesystem::path bench_dir, T& method)
-{
-    traverse_recursive(parent, "", parent, bench_dir, method);
-}
-
-void replaceIncludes(
-    [[maybe_unused]] boost::property_tree::ptree const& parent,
-    [[maybe_unused]] boost::property_tree::ptree::path_type const& child_path,
-    boost::property_tree::ptree& child,
-    std::filesystem::path const& bench_dir)
-{
-    using boost::property_tree::ptree;
-    for (auto& [key, tree] : child)
-    {
-        if (key == "include")
-        {
-            auto filename = tree.get<std::string>("<xmlattr>.file");
-            if (auto const filepath = std::filesystem::path(filename);
-                filepath.is_relative())
-            {
-                filename = (bench_dir / filepath).string();
-            }
-            INFO("Including {:s} into project file.", filename);
-
-            ptree include_tree;
-            read_xml(filename, include_tree,
-                     boost::property_tree::xml_parser::no_comments |
-                         boost::property_tree::xml_parser::trim_whitespace);
-
-            // Can only insert subtree at child
-            auto& tmp_tree = child.put_child("include", include_tree);
-
-            // Move subtree above child
-            std::move(tmp_tree.begin(), tmp_tree.end(), back_inserter(child));
-
-            // Erase child
-            child.erase("include");
-
-            // There can only be one include under a parent element!
-            break;
-        }
-    }
-}
-
-// Applies a patch file to the prj content in prj_stream.
-void patchStream(std::string patch_file, std::stringstream& prj_stream)
-{
-    auto patch = xmlParseFile(patch_file.c_str());
-    if (patch == NULL)
-    {
-        OGS_FATAL("Error reading XML diff file {:s}.", patch_file);
-    }
-
-    auto doc =
-        xmlParseMemory(prj_stream.str().c_str(), prj_stream.str().size());
-    if (doc == NULL)
-    {
-        OGS_FATAL("Error reading project file from memory.");
-    }
-
-    auto node = xmlDocGetRootElement(patch);
-    int rc = 0;
-    for (node = node ? node->children : NULL; node; node = node->next)
-    {
-        if (node->type != XML_ELEMENT_NODE)
-        {
-            continue;
-        }
-
-        if (!strcmp(reinterpret_cast<const char*>(node->name), "add"))
-        {
-            rc = xml_patch_add(doc, node);
-        }
-        else if (!strcmp(reinterpret_cast<const char*>(node->name), "replace"))
-        {
-            rc = xml_patch_replace(doc, node);
-        }
-        else if (!strcmp(reinterpret_cast<const char*>(node->name), "remove"))
-        {
-            rc = xml_patch_remove(doc, node);
-        }
-        else
-        {
-            OGS_FATAL(
-                "Error while patching prj file with patch file {:}. Only "
-                "'add', 'replace' and 'remove' elements are allowed! Got an "
-                "element '{:s}' on line {:d}.",
-                patch_file, node->name, node->line);
-        }
-
-        if (rc)
-        {
-            OGS_FATAL(
-                "Error while patching prj file with patch file {:}. Error in "
-                "element '{:s}' on line {:d}.",
-                patch_file, node->name, node->line);
-        }
-    }
-
-    xmlChar* xmlbuff;
-    int buffersize;
-    xmlDocDumpMemory(doc, &xmlbuff, &buffersize);
-    prj_stream.str("");  // Empty stream
-    prj_stream << xmlbuff;
-
-    xmlFree(xmlbuff);
-    xmlFreeDoc(doc);
-    xmlFreeDoc(patch);
-}
-
-// Will set prj_file to the actual .prj file and returns the final prj file
-// content in prj_stream.
-void readAndPatchPrj(std::stringstream& prj_stream, std::string& prj_file,
-                     std::vector<std::string> patch_files)
-{
-    // Extract base project file path if an xml (patch) file is given as prj
-    // file and it contains the base_file attribute.
-    if (BaseLib::getFileExtension(prj_file) == ".xml")
-    {
-        if (!patch_files.empty())
-        {
-            OGS_FATAL(
-                "It is not allowed to specify additional patch files "
-                "if a patch file was already specified as the "
-                "prj-file.");
-        }
-        auto patch = xmlParseFile(prj_file.c_str());
-        auto node = xmlDocGetRootElement(patch);
-        auto base_file = xmlGetProp(node, (const xmlChar*)"base_file");
-        if (base_file == nullptr)
-        {
-            OGS_FATAL(
-                "Error reading base prj file (base_file attribute) in given "
-                "patch file {:s}.",
-                prj_file);
-        }
-        patch_files = {prj_file};
-        std::stringstream ss;
-        ss << base_file;
-        prj_file = BaseLib::joinPaths(BaseLib::extractPath(prj_file), ss.str());
-    }
-
-    // read base prj file into stream
-    if (std::ifstream file(prj_file); file)
-    {
-        prj_stream << file.rdbuf();
-    }
-    else
-    {
-        if (!BaseLib::IsFileExisting(prj_file))
-        {
-            ERR("File {:s} does not exist.", prj_file);
-        }
-        DBUG("Stream state flags: {:b}.", file.rdstate());
-        OGS_FATAL("Could not open project file '{:s}' for reading.", prj_file);
-    }
-
-    // apply xml patches to stream
-    if (!patch_files.empty())
-    {
-        for (const auto& patch_file : patch_files)
-        {
-            patchStream(patch_file, prj_stream);
-        }
-        xmlCleanupParser();
-    }
-}
-
-ConfigTree makeConfigTree(const std::string& filepath,
-                          const bool be_ruthless,
+ConfigTree makeConfigTree(const std::string& filepath, const bool be_ruthless,
                           const std::string& toplevel_tag,
-                          const std::vector<std::string>& patch_files)
+                          std::stringstream& prj_stream)
 {
-    std::string prj_file = filepath;
-    std::stringstream prj_stream;
-    readAndPatchPrj(prj_stream, prj_file, patch_files);
-
     ConfigTree::PTree ptree;
 
     // note: Trimming whitespace and ignoring comments is crucial in order
@@ -233,12 +31,6 @@ ConfigTree makeConfigTree(const std::string& filepath,
         read_xml(prj_stream, ptree,
                  boost::property_tree::xml_parser::no_comments |
                      boost::property_tree::xml_parser::trim_whitespace);
-
-        if (toplevel_tag == "OpenGeoSysProject")
-        {
-            traverse(ptree, std::filesystem::path(prj_file).parent_path(),
-                     replaceIncludes);
-        }
     }
     catch (boost::property_tree::xml_parser_error const& e)
     {
@@ -259,4 +51,22 @@ ConfigTree makeConfigTree(const std::string& filepath,
               filepath);
 }
 
+ConfigTree makeConfigTreeFromFile(const std::string& filepath,
+                                  const bool be_ruthless,
+                                  const std::string& toplevel_tag)
+{
+    std::ifstream file(filepath);
+    if (file)
+    {
+        std::stringstream buffer;
+        buffer << file.rdbuf();
+
+        return makeConfigTree(filepath, be_ruthless, toplevel_tag, buffer);
+    }
+    else
+    {
+        OGS_FATAL("Could not read from file {:s}!", filepath);
+    }
+}
+
 }  // namespace BaseLib
diff --git a/BaseLib/ConfigTreeUtil.h b/BaseLib/ConfigTreeUtil.h
index 08d48923e18b84fcfe5b700fd1458159217d121a..1026272d2229e35671df53f41c7fb09d87a4e8d5 100644
--- a/BaseLib/ConfigTreeUtil.h
+++ b/BaseLib/ConfigTreeUtil.h
@@ -11,7 +11,6 @@
 #pragma once
 
 #include "ConfigTree.h"
-#include "Logging.h"
 
 namespace BaseLib
 {
@@ -21,7 +20,7 @@ namespace BaseLib
  * \param be_ruthless  see ConfigTreeTopLevel::ConfigTreeTopLevel()
  * \param toplevel_tag name of the outermost tag in the XML file. The returned
  *                     ConfigTree is rooted one level below that tag.
- * \param patch_files  optional vector of strings with patch file paths.
+ * \param prj_stream   content of the (pre-processed) prj file.
  *
  * The parameter \c toplevel_tag is provided for compatibility with our existing
  * configuration files whose toplevel tags are written in camel case, which
@@ -45,6 +44,10 @@ namespace BaseLib
 ConfigTree makeConfigTree(std::string const& filepath,
                           bool const be_ruthless,
                           std::string const& toplevel_tag,
-                          const std::vector<std::string>& patch_files = {});
+                          std::stringstream& prj_stream);
+
+ConfigTree makeConfigTreeFromFile(const std::string& filepath,
+                                  const bool be_ruthless,
+                                  const std::string& toplevel_tag);
 
 }  // namespace BaseLib
diff --git a/BaseLib/PrjProcessing.cpp b/BaseLib/PrjProcessing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..84702330f4594c7e7d4bb67ae7a7e95983ba4fa0
--- /dev/null
+++ b/BaseLib/PrjProcessing.cpp
@@ -0,0 +1,310 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 "PrjProcessing.h"
+
+#include <libxml/parser.h>
+#include <libxml/tree.h>
+#include <xml_patch.h>
+
+#include <filesystem>
+#include <fstream>
+
+#include "Error.h"
+#include "FileTools.h"
+#include "Logging.h"
+
+namespace BaseLib
+{
+void traverseIncludes(xmlDoc* doc, xmlNode* node,
+                      std::filesystem::path const& prj_dir)
+{
+    xmlNode* cur_node = nullptr;
+    for (cur_node = node; cur_node; cur_node = cur_node->next)
+    {
+        if (cur_node->type != XML_ELEMENT_NODE)
+        {
+            continue;
+        }
+        if (xmlStrEqual(cur_node->name, xmlCharStrdup("include")))
+        {
+            auto include_file_char_pointer =
+                xmlGetProp(cur_node, xmlCharStrdup("file"));
+            if (include_file_char_pointer == nullptr)
+            {
+                OGS_FATAL(
+                    "Error while processing includes in prj file. Error in "
+                    "element '{:s}' on line {:d}: no file attribute given!",
+                    cur_node->name, cur_node->line);
+            }
+            auto filename_length = xmlStrlen(include_file_char_pointer);
+            std::string filename(
+                reinterpret_cast<char*>(include_file_char_pointer),
+                filename_length);
+            if (auto const filepath = std::filesystem::path(filename);
+                filepath.is_relative())
+            {
+                filename = (prj_dir / filepath).string();
+            }
+
+            if (!std::filesystem::exists(filename))
+            {
+                OGS_FATAL(
+                    "Error while processing includes in prj file. Error in "
+                    "element '{:s}' on line {:d}: Include file is not "
+                    "existing: "
+                    "{:s}!",
+                    cur_node->name, cur_node->line, include_file_char_pointer);
+            }
+            INFO("Including {:s} into project file.", filename);
+
+            const std::ifstream input_stream(filename, std::ios_base::binary);
+            if (input_stream.fail())
+            {
+                OGS_FATAL("Failed to open file {s}!", filename);
+            }
+            std::stringstream buffer;
+            buffer << input_stream.rdbuf();
+            const std::string xml = buffer.str();
+
+            xmlNodePtr pNewNode = nullptr;
+            xmlParseInNodeContext(cur_node->parent, xml.c_str(),
+                                  (int)xml.length(), 0, &pNewNode);
+            if (pNewNode != nullptr)
+            {
+                // add new xml node to parent
+                xmlNode* pChild = pNewNode;
+                while (pChild != nullptr)
+                {
+                    xmlAddChild(cur_node->parent, xmlCopyNode(pChild, 1));
+                    pChild = pChild->next;
+                }
+                xmlFreeNode(pNewNode);
+            }
+
+            // Cleanup and continue on next node
+            auto next_node = cur_node->next;
+            xmlUnlinkNode(cur_node);
+            xmlFreeNode(cur_node);
+            cur_node = next_node;
+        }
+        traverseIncludes(doc, cur_node->children, prj_dir);
+    }
+    xmlFreeNode(cur_node);
+}
+
+void replaceIncludes(std::stringstream& prj_stream,
+                     std::filesystem::path const& prj_dir)
+{
+    auto doc =
+        xmlParseMemory(prj_stream.str().c_str(), prj_stream.str().size());
+    if (doc == nullptr)
+    {
+        OGS_FATAL("Error reading project file from memory.");
+    }
+
+    auto root_node = xmlDocGetRootElement(doc);
+    traverseIncludes(doc, root_node, prj_dir);
+
+    xmlChar* xmlbuff;
+    int buffersize;
+    xmlDocDumpMemory(doc, &xmlbuff, &buffersize);
+    prj_stream.str("");  // Empty stream
+    prj_stream << xmlbuff;
+
+    xmlFree(xmlbuff);
+    xmlFreeDoc(doc);
+}
+
+// Applies a patch file to the prj content in prj_stream.
+void patchStream(std::string const& patch_file, std::stringstream& prj_stream,
+                 bool after_includes = false)
+{
+    auto patch = xmlParseFile(patch_file.c_str());
+    if (patch == nullptr)
+    {
+        OGS_FATAL("Error reading XML diff file {:s}.", patch_file);
+    }
+
+    auto doc =
+        xmlParseMemory(prj_stream.str().c_str(), prj_stream.str().size());
+    if (doc == nullptr)
+    {
+        OGS_FATAL("Error reading project file from memory.");
+    }
+
+    auto node = xmlDocGetRootElement(patch);
+    int rc = 0;
+    for (node = node ? node->children : nullptr; node; node = node->next)
+    {
+        if (node->type != XML_ELEMENT_NODE)
+        {
+            continue;
+        }
+        bool node_after_includes = false;
+        xmlAttr* attribute = node->properties;
+        while (attribute)
+        {
+            // Check for after_includes-attribute
+            xmlChar* value =
+                xmlNodeListGetString(node->doc, attribute->children, 1);
+            if (xmlStrEqual(attribute->name, xmlCharStrdup("after_includes")) &&
+                xmlStrEqual(value, xmlCharStrdup("true")))
+            {
+                node_after_includes = true;
+            }
+
+            xmlFree(value);
+            attribute = attribute->next;
+        }
+
+        if (after_includes != node_after_includes)
+        {
+            continue;
+        }
+
+        if (xmlStrEqual(node->name, xmlCharStrdup("add")))
+        {
+            rc = xml_patch_add(doc, node);
+        }
+        else if (xmlStrEqual(node->name, xmlCharStrdup("replace")))
+        {
+            rc = xml_patch_replace(doc, node);
+        }
+        else if (xmlStrEqual(node->name, xmlCharStrdup("remove")))
+        {
+            rc = xml_patch_remove(doc, node);
+        }
+        else
+        {
+            OGS_FATAL(
+                "Error while patching prj file with patch file {:}. Only "
+                "'add', 'replace' and 'remove' elements are allowed! Got an "
+                "element '{:s}' on line {:d}.",
+                patch_file, node->name, node->line);
+        }
+
+        if (rc)
+        {
+            OGS_FATAL(
+                "Error while patching prj file with patch file {:}. Error in "
+                "element '{:s}' on line {:d}.",
+                patch_file, node->name, node->line);
+        }
+    }
+
+    xmlChar* xmlbuff;
+    int buffersize;
+    xmlDocDumpMemory(doc, &xmlbuff, &buffersize);
+    prj_stream.str("");  // Empty stream
+    prj_stream << xmlbuff;
+
+    xmlFree(xmlbuff);
+    xmlFreeDoc(doc);
+    xmlFreeDoc(patch);
+}
+
+// Will set prj_file to the actual .prj file and returns the final prj file
+// content in prj_stream.
+void readAndPatchPrj(std::stringstream& prj_stream, std::string& prj_file,
+                     std::vector<std::string>& patch_files)
+{
+    // Extract base project file path if an xml (patch) file is given as prj
+    // file and it contains the base_file attribute.
+    if (BaseLib::getFileExtension(prj_file) == ".xml")
+    {
+        if (!patch_files.empty())
+        {
+            OGS_FATAL(
+                "It is not allowed to specify additional patch files "
+                "if a patch file was already specified as the "
+                "prj-file.");
+        }
+        auto patch = xmlParseFile(prj_file.c_str());
+        auto node = xmlDocGetRootElement(patch);
+        auto base_file = xmlGetProp(node, (const xmlChar*)"base_file");
+        if (base_file == nullptr)
+        {
+            OGS_FATAL(
+                "Error reading base prj file (base_file attribute) in given "
+                "patch file {:s}.",
+                prj_file);
+        }
+        patch_files = {prj_file};
+        std::stringstream ss;
+        ss << base_file;
+        prj_file = BaseLib::joinPaths(BaseLib::extractPath(prj_file), ss.str());
+    }
+
+    // read base prj file into stream
+    if (std::ifstream file(prj_file); file)
+    {
+        prj_stream << file.rdbuf();
+    }
+    else
+    {
+        if (!BaseLib::IsFileExisting(prj_file))
+        {
+            ERR("File {:s} does not exist.", prj_file);
+        }
+        DBUG("Stream state flags: {:b}.", file.rdstate());
+        OGS_FATAL("Could not open project file '{:s}' for reading.", prj_file);
+    }
+
+    // apply xml patches to stream
+    for (const auto& patch_file : patch_files)
+    {
+        patchStream(patch_file, prj_stream);
+    }
+}
+
+void prepareProjectFile(std::stringstream& prj_stream,
+                        const std::string& filepath,
+                        const std::vector<std::string>& patch_files,
+                        bool write_prj, const std::string& out_directory)
+{
+    std::string prj_file = filepath;
+
+    std::vector<std::string> patch_files_copy = patch_files;
+    readAndPatchPrj(prj_stream, prj_file, patch_files_copy);
+    // LB TODO later: replace canonical with absolute when a mesh input dir
+    // ogs parameter is implemented.
+    replaceIncludes(prj_stream,
+                    std::filesystem::canonical(std::filesystem::path(prj_file))
+                        .parent_path());
+    // re-apply xml patches to stream
+    for (const auto& patch_file : patch_files_copy)
+    {
+        patchStream(patch_file, prj_stream, true);
+    }
+
+    if (write_prj)
+    {
+        // pretty-print
+        xmlKeepBlanksDefault(0);
+
+        // The following two lines should set indenting to 4 spaces but it does
+        // not work. 2 spaces are default.
+        //
+        // xmlThrDefIndentTreeOutput(1);
+        // xmlThrDefTreeIndentString("    "); // 4 spaces indent
+        auto doc =
+            xmlParseMemory(prj_stream.str().c_str(), prj_stream.str().size());
+        auto prj_out = (std::filesystem::path(out_directory) /
+                        std::filesystem::path(filepath).stem())
+                           .string() +
+                       "_processed.prj";
+        xmlSaveFormatFileEnc(prj_out.c_str(), doc, "utf-8", 1);
+        INFO("Processed project file written to {:s}.", prj_out);
+        xmlFreeDoc(doc);
+    }
+    xmlCleanupParser();
+}
+}  // namespace BaseLib
diff --git a/BaseLib/PrjProcessing.h b/BaseLib/PrjProcessing.h
new file mode 100644
index 0000000000000000000000000000000000000000..138ee7d0faa5f7acc6db61c22fec98362731e86d
--- /dev/null
+++ b/BaseLib/PrjProcessing.h
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, 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 <sstream>
+#include <vector>
+
+namespace BaseLib
+{
+/**
+ * \brief Applies includes and patch files to project file.
+ *
+ * \param prj_stream    The processed prj as a stringstream.
+ * \param filepath      The prj file.
+ * \param patch_files   Optional patch files.
+ * \param write_prj     Write the processed project file to disk?
+ * \param out_directory The directory to write the processed file to.
+ */
+void prepareProjectFile(std::stringstream& prj_stream,
+                        const std::string& filepath,
+                        const std::vector<std::string>& patch_files,
+                        bool write_prj, const std::string& out_directory);
+}  // namespace BaseLib
diff --git a/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp b/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp
index 9ec613610487604b6da5f8dcecc19de6d92d8178..a3de9441278f7190d341e430787c9de0541b8b81 100644
--- a/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp
+++ b/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp
@@ -38,7 +38,7 @@ BoostXmlGmlInterface::BoostXmlGmlInterface(GeoLib::GEOObjects& geo_objs)
 bool BoostXmlGmlInterface::readFile(const std::string& fname)
 {
     //! \todo Reading geometries is always strict.
-    auto doc = BaseLib::makeConfigTree(fname, true, "OpenGeoSysGLI");
+    auto doc = BaseLib::makeConfigTreeFromFile(fname, true, "OpenGeoSysGLI");
 
     // ignore attributes related to XML schema
     doc.ignoreConfigAttribute("xmlns:xsi");
diff --git a/ProcessLib/SteadyStateDiffusion/Tests.cmake b/ProcessLib/SteadyStateDiffusion/Tests.cmake
index c7c7c6de2f6fe13871f7db467340b8057f1904c1..3a8db7a2224f3af20bd9df5782d9161dfdd02974 100644
--- a/ProcessLib/SteadyStateDiffusion/Tests.cmake
+++ b/ProcessLib/SteadyStateDiffusion/Tests.cmake
@@ -4,13 +4,21 @@ foreach(mesh_size 1e0 1e1 1e2 1e3)
         NAME SteadyStateDiffusion_cube_1x1x1_${mesh_size}
         PATH Elliptic/cube_1x1x1_SteadyStateDiffusion
         EXECUTABLE ogs
-        EXECUTABLE_ARGS cube_${mesh_size}.xml
+        EXECUTABLE_ARGS cube_${mesh_size}.xml --write-prj
         TESTER vtkdiff
         REQUIREMENTS NOT OGS_USE_MPI
         DIFF_DATA
         cube_1x1x1_hex_${mesh_size}.vtu cube_${mesh_size}_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure 1e-15 1e-15
     )
 
+    if(TEST SteadyStateDiffusion_cube_1x1x1_${mesh_size} AND DIFF_TOOL_PATH)
+        set(_processed_path Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_${mesh_size}_processed.prj)
+        add_test(NAME SteadyStateDiffusion_cube_1x1x1_${mesh_size}_prj_diff
+            COMMAND ${DIFF_TOOL_PATH} ${Data_SOURCE_DIR}/${_processed_path} ${Data_BINARY_DIR}/${_processed_path})
+        set_tests_properties(SteadyStateDiffusion_cube_1x1x1_${mesh_size}_prj_diff
+            PROPERTIES LABELS "default" DEPENDS SteadyStateDiffusion_cube_1x1x1_${mesh_size})
+    endif()
+
     AddTest(
         NAME SteadyStateDiffusion_cube_1x1x1_${mesh_size}_Newton
         PATH Elliptic/cube_1x1x1_SteadyStateDiffusion
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0.prj b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0.prj
index ac38a2b8266c89fda3be2854d440943bf28b6b12..cebf20775181a10c671c40a1ac474b46b02f37c2 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0.prj
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0.prj
@@ -38,15 +38,7 @@
                 </time_stepping>
             </process>
         </processes>
-        <output>
-            <type>VTK</type>
-            <prefix>cube_1e0</prefix>
-            <variables>
-                <variable> pressure </variable>
-                <variable> v      </variable>
-            </variables>
-            <suffix>_ts_{:timestep}_t_{:time}</suffix>
-        </output>
+        <include file="cube_output.xml" />
     </time_loop>
     <parameters>
         <parameter>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_processed.prj b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_processed.prj
new file mode 100644
index 0000000000000000000000000000000000000000..3a18d6e916657aa3a54528c43ecfe701446d34ee
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_processed.prj
@@ -0,0 +1,129 @@
+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysProject>
+  <mesh>cube_1x1x1_hex_1e0.vtu</mesh>
+  <geometry>cube_1x1x1.gml</geometry>
+  <processes>
+    <process>
+      <name>SteadyStateDiffusion</name>
+      <type>STEADY_STATE_DIFFUSION</type>
+      <integration_order>2</integration_order>
+      <process_variables>
+        <process_variable>pressure</process_variable>
+      </process_variables>
+      <secondary_variables>
+        <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+      </secondary_variables>
+    </process>
+  </processes>
+  <media>
+    <medium id="0">
+      <properties>
+        <property>
+          <name>diffusion</name>
+          <type>Constant</type>
+          <value>1.0</value>
+        </property>
+        <property>
+          <name>reference_temperature</name>
+          <type>Constant</type>
+          <value>293.15</value>
+        </property>
+      </properties>
+    </medium>
+  </media>
+  <time_loop>
+    <processes>
+      <process ref="SteadyStateDiffusion">
+        <nonlinear_solver>basic_picard</nonlinear_solver>
+        <convergence_criterion>
+          <type>DeltaX</type>
+          <norm_type>NORM2</norm_type>
+          <abstol>1.e-6</abstol>
+        </convergence_criterion>
+        <time_discretization>
+          <type>BackwardEuler</type>
+        </time_discretization>
+        <time_stepping>
+          <type>SingleStep</type>
+        </time_stepping>
+      </process>
+    </processes>
+    <output>
+      <type>VTK</type>
+      <prefix>cube_1e0</prefix>
+      <variables>
+        <variable> pressure </variable>
+        <variable> v      </variable>
+      </variables>
+      <suffix>_ts_{:timestep}_t_{:time}</suffix>
+    </output>
+  </time_loop>
+  <parameters>
+    <parameter>
+      <name>K</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p0</name>
+      <type>Constant</type>
+      <value>0</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_left</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_right</name>
+      <type>Constant</type>
+      <value>-1</value>
+    </parameter>
+  </parameters>
+  <process_variables>
+    <process_variable>
+      <name>pressure</name>
+      <components>1</components>
+      <order>1</order>
+      <initial_condition>p0</initial_condition>
+      <boundary_conditions>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>left</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_left</parameter>
+        </boundary_condition>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>right</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_right</parameter>
+        </boundary_condition>
+      </boundary_conditions>
+    </process_variable>
+  </process_variables>
+  <nonlinear_solvers>
+    <nonlinear_solver>
+      <name>basic_picard</name>
+      <type>Picard</type>
+      <max_iter>10</max_iter>
+      <linear_solver>general_linear_solver</linear_solver>
+    </nonlinear_solver>
+  </nonlinear_solvers>
+  <linear_solvers>
+    <linear_solver>
+      <name>general_linear_solver</name>
+      <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+      <eigen>
+        <solver_type>CG</solver_type>
+        <precon_type>DIAGONAL</precon_type>
+        <max_iteration_step>10000</max_iteration_step>
+        <error_tolerance>1e-16</error_tolerance>
+      </eigen>
+      <petsc>
+        <prefix>gw</prefix>
+        <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+      </petsc>
+    </linear_solver>
+  </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1.xml b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1.xml
index ae07b72e625af23372c10c0316b5568605f11cc2..0a48e164967723c5a2e3efbaedbbda8482a24bed 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1.xml
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1.xml
@@ -1,5 +1,5 @@
 <?xml version="1.0" encoding="ISO-8859-1"?>
 <OpenGeoSysProjectDiff base_file="cube_1e0.prj">
     <replace sel="/*/mesh/text()">cube_1x1x1_hex_1e1.vtu</replace>
-    <replace sel="/*/time_loop/output/prefix/text()">cube_1e1</replace>
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">cube_1e1</replace>
 </OpenGeoSysProjectDiff>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1_processed.prj b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1_processed.prj
new file mode 100644
index 0000000000000000000000000000000000000000..3ecce15d6b27b8a687901cd1db0626640996a626
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e1_processed.prj
@@ -0,0 +1,129 @@
+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysProject>
+  <mesh>cube_1x1x1_hex_1e1.vtu</mesh>
+  <geometry>cube_1x1x1.gml</geometry>
+  <processes>
+    <process>
+      <name>SteadyStateDiffusion</name>
+      <type>STEADY_STATE_DIFFUSION</type>
+      <integration_order>2</integration_order>
+      <process_variables>
+        <process_variable>pressure</process_variable>
+      </process_variables>
+      <secondary_variables>
+        <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+      </secondary_variables>
+    </process>
+  </processes>
+  <media>
+    <medium id="0">
+      <properties>
+        <property>
+          <name>diffusion</name>
+          <type>Constant</type>
+          <value>1.0</value>
+        </property>
+        <property>
+          <name>reference_temperature</name>
+          <type>Constant</type>
+          <value>293.15</value>
+        </property>
+      </properties>
+    </medium>
+  </media>
+  <time_loop>
+    <processes>
+      <process ref="SteadyStateDiffusion">
+        <nonlinear_solver>basic_picard</nonlinear_solver>
+        <convergence_criterion>
+          <type>DeltaX</type>
+          <norm_type>NORM2</norm_type>
+          <abstol>1.e-6</abstol>
+        </convergence_criterion>
+        <time_discretization>
+          <type>BackwardEuler</type>
+        </time_discretization>
+        <time_stepping>
+          <type>SingleStep</type>
+        </time_stepping>
+      </process>
+    </processes>
+    <output>
+      <type>VTK</type>
+      <prefix>cube_1e1</prefix>
+      <variables>
+        <variable> pressure </variable>
+        <variable> v      </variable>
+      </variables>
+      <suffix>_ts_{:timestep}_t_{:time}</suffix>
+    </output>
+  </time_loop>
+  <parameters>
+    <parameter>
+      <name>K</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p0</name>
+      <type>Constant</type>
+      <value>0</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_left</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_right</name>
+      <type>Constant</type>
+      <value>-1</value>
+    </parameter>
+  </parameters>
+  <process_variables>
+    <process_variable>
+      <name>pressure</name>
+      <components>1</components>
+      <order>1</order>
+      <initial_condition>p0</initial_condition>
+      <boundary_conditions>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>left</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_left</parameter>
+        </boundary_condition>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>right</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_right</parameter>
+        </boundary_condition>
+      </boundary_conditions>
+    </process_variable>
+  </process_variables>
+  <nonlinear_solvers>
+    <nonlinear_solver>
+      <name>basic_picard</name>
+      <type>Picard</type>
+      <max_iter>10</max_iter>
+      <linear_solver>general_linear_solver</linear_solver>
+    </nonlinear_solver>
+  </nonlinear_solvers>
+  <linear_solvers>
+    <linear_solver>
+      <name>general_linear_solver</name>
+      <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+      <eigen>
+        <solver_type>CG</solver_type>
+        <precon_type>DIAGONAL</precon_type>
+        <max_iteration_step>10000</max_iteration_step>
+        <error_tolerance>1e-16</error_tolerance>
+      </eigen>
+      <petsc>
+        <prefix>gw</prefix>
+        <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+      </petsc>
+    </linear_solver>
+  </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2.xml b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2.xml
index c371070e2780402a1d29d6ad71feb4c8dbbd56b5..d4adff5c134a576750e8567a42a6b1061b58a432 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2.xml
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2.xml
@@ -1,5 +1,5 @@
 <?xml version="1.0" encoding="ISO-8859-1"?>
 <OpenGeoSysProjectDiff base_file="cube_1e0.prj">
     <replace sel="/*/mesh/text()">cube_1x1x1_hex_1e2.vtu</replace>
-    <replace sel="/*/time_loop/output/prefix/text()">cube_1e2</replace>
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">cube_1e2</replace>
 </OpenGeoSysProjectDiff>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2_processed.prj b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2_processed.prj
new file mode 100644
index 0000000000000000000000000000000000000000..bb1f7787253c76cdf64cd2a0a21d35f418c6852c
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e2_processed.prj
@@ -0,0 +1,129 @@
+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysProject>
+  <mesh>cube_1x1x1_hex_1e2.vtu</mesh>
+  <geometry>cube_1x1x1.gml</geometry>
+  <processes>
+    <process>
+      <name>SteadyStateDiffusion</name>
+      <type>STEADY_STATE_DIFFUSION</type>
+      <integration_order>2</integration_order>
+      <process_variables>
+        <process_variable>pressure</process_variable>
+      </process_variables>
+      <secondary_variables>
+        <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+      </secondary_variables>
+    </process>
+  </processes>
+  <media>
+    <medium id="0">
+      <properties>
+        <property>
+          <name>diffusion</name>
+          <type>Constant</type>
+          <value>1.0</value>
+        </property>
+        <property>
+          <name>reference_temperature</name>
+          <type>Constant</type>
+          <value>293.15</value>
+        </property>
+      </properties>
+    </medium>
+  </media>
+  <time_loop>
+    <processes>
+      <process ref="SteadyStateDiffusion">
+        <nonlinear_solver>basic_picard</nonlinear_solver>
+        <convergence_criterion>
+          <type>DeltaX</type>
+          <norm_type>NORM2</norm_type>
+          <abstol>1.e-6</abstol>
+        </convergence_criterion>
+        <time_discretization>
+          <type>BackwardEuler</type>
+        </time_discretization>
+        <time_stepping>
+          <type>SingleStep</type>
+        </time_stepping>
+      </process>
+    </processes>
+    <output>
+      <type>VTK</type>
+      <prefix>cube_1e2</prefix>
+      <variables>
+        <variable> pressure </variable>
+        <variable> v      </variable>
+      </variables>
+      <suffix>_ts_{:timestep}_t_{:time}</suffix>
+    </output>
+  </time_loop>
+  <parameters>
+    <parameter>
+      <name>K</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p0</name>
+      <type>Constant</type>
+      <value>0</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_left</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_right</name>
+      <type>Constant</type>
+      <value>-1</value>
+    </parameter>
+  </parameters>
+  <process_variables>
+    <process_variable>
+      <name>pressure</name>
+      <components>1</components>
+      <order>1</order>
+      <initial_condition>p0</initial_condition>
+      <boundary_conditions>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>left</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_left</parameter>
+        </boundary_condition>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>right</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_right</parameter>
+        </boundary_condition>
+      </boundary_conditions>
+    </process_variable>
+  </process_variables>
+  <nonlinear_solvers>
+    <nonlinear_solver>
+      <name>basic_picard</name>
+      <type>Picard</type>
+      <max_iter>10</max_iter>
+      <linear_solver>general_linear_solver</linear_solver>
+    </nonlinear_solver>
+  </nonlinear_solvers>
+  <linear_solvers>
+    <linear_solver>
+      <name>general_linear_solver</name>
+      <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+      <eigen>
+        <solver_type>CG</solver_type>
+        <precon_type>DIAGONAL</precon_type>
+        <max_iteration_step>10000</max_iteration_step>
+        <error_tolerance>1e-16</error_tolerance>
+      </eigen>
+      <petsc>
+        <prefix>gw</prefix>
+        <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+      </petsc>
+    </linear_solver>
+  </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3.xml b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3.xml
index daf9b37ee591ee288dbf0d2aab71246547c00b68..2d6c62af6c16f3c334a7b17b8ebd30e9b368d687 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3.xml
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3.xml
@@ -1,5 +1,5 @@
 <?xml version="1.0" encoding="ISO-8859-1"?>
 <OpenGeoSysProjectDiff base_file="cube_1e0.prj">
     <replace sel="/*/mesh/text()">cube_1x1x1_hex_1e3.vtu</replace>
-    <replace sel="/*/time_loop/output/prefix/text()">cube_1e3</replace>
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">cube_1e3</replace>
 </OpenGeoSysProjectDiff>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3_processed.prj b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3_processed.prj
new file mode 100644
index 0000000000000000000000000000000000000000..da41f5a1bd2b4ce353e6b1fbf40187b696c91b57
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e3_processed.prj
@@ -0,0 +1,129 @@
+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysProject>
+  <mesh>cube_1x1x1_hex_1e3.vtu</mesh>
+  <geometry>cube_1x1x1.gml</geometry>
+  <processes>
+    <process>
+      <name>SteadyStateDiffusion</name>
+      <type>STEADY_STATE_DIFFUSION</type>
+      <integration_order>2</integration_order>
+      <process_variables>
+        <process_variable>pressure</process_variable>
+      </process_variables>
+      <secondary_variables>
+        <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+      </secondary_variables>
+    </process>
+  </processes>
+  <media>
+    <medium id="0">
+      <properties>
+        <property>
+          <name>diffusion</name>
+          <type>Constant</type>
+          <value>1.0</value>
+        </property>
+        <property>
+          <name>reference_temperature</name>
+          <type>Constant</type>
+          <value>293.15</value>
+        </property>
+      </properties>
+    </medium>
+  </media>
+  <time_loop>
+    <processes>
+      <process ref="SteadyStateDiffusion">
+        <nonlinear_solver>basic_picard</nonlinear_solver>
+        <convergence_criterion>
+          <type>DeltaX</type>
+          <norm_type>NORM2</norm_type>
+          <abstol>1.e-6</abstol>
+        </convergence_criterion>
+        <time_discretization>
+          <type>BackwardEuler</type>
+        </time_discretization>
+        <time_stepping>
+          <type>SingleStep</type>
+        </time_stepping>
+      </process>
+    </processes>
+    <output>
+      <type>VTK</type>
+      <prefix>cube_1e3</prefix>
+      <variables>
+        <variable> pressure </variable>
+        <variable> v      </variable>
+      </variables>
+      <suffix>_ts_{:timestep}_t_{:time}</suffix>
+    </output>
+  </time_loop>
+  <parameters>
+    <parameter>
+      <name>K</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p0</name>
+      <type>Constant</type>
+      <value>0</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_left</name>
+      <type>Constant</type>
+      <value>1</value>
+    </parameter>
+    <parameter>
+      <name>p_Dirichlet_right</name>
+      <type>Constant</type>
+      <value>-1</value>
+    </parameter>
+  </parameters>
+  <process_variables>
+    <process_variable>
+      <name>pressure</name>
+      <components>1</components>
+      <order>1</order>
+      <initial_condition>p0</initial_condition>
+      <boundary_conditions>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>left</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_left</parameter>
+        </boundary_condition>
+        <boundary_condition>
+          <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+          <geometry>right</geometry>
+          <type>Dirichlet</type>
+          <parameter>p_Dirichlet_right</parameter>
+        </boundary_condition>
+      </boundary_conditions>
+    </process_variable>
+  </process_variables>
+  <nonlinear_solvers>
+    <nonlinear_solver>
+      <name>basic_picard</name>
+      <type>Picard</type>
+      <max_iter>10</max_iter>
+      <linear_solver>general_linear_solver</linear_solver>
+    </nonlinear_solver>
+  </nonlinear_solvers>
+  <linear_solvers>
+    <linear_solver>
+      <name>general_linear_solver</name>
+      <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+      <eigen>
+        <solver_type>CG</solver_type>
+        <precon_type>DIAGONAL</precon_type>
+        <max_iteration_step>10000</max_iteration_step>
+        <error_tolerance>1e-16</error_tolerance>
+      </eigen>
+      <petsc>
+        <prefix>gw</prefix>
+        <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+      </petsc>
+    </linear_solver>
+  </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_output.xml b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_output.xml
new file mode 100644
index 0000000000000000000000000000000000000000..1f2ec21c9099073ad5d964c29c4190443c1e5893
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_output.xml
@@ -0,0 +1,9 @@
+<output>
+    <type>VTK</type>
+    <prefix>cube_1e0</prefix>
+    <variables>
+        <variable> pressure </variable>
+        <variable> v      </variable>
+    </variables>
+    <suffix>_ts_{:timestep}_t_{:time}</suffix>
+</output>
diff --git a/web/content/docs/userguide/basics/cli-arguments.md b/web/content/docs/userguide/basics/cli-arguments.md
index d03e3bee88cbd70f63a11eff40c8b9d037674efc..8b72094d855fb2d9cf4d92a5c08a332132526115 100644
--- a/web/content/docs/userguide/basics/cli-arguments.md
+++ b/web/content/docs/userguide/basics/cli-arguments.md
@@ -30,6 +30,10 @@ $ ogs --help
    -o <PATH>,  --output-directory <PATH>
      the output directory to write to
 
+   --write-prj
+     Writes processed project file to output path /
+     [prj_base_name]_processed.prj.
+
    -p <>,  --xml-patch <>  (accepted multiple times)
      the xml patch file(s) which is (are) applied (in the given order) to
      the PROJECT_FILE
diff --git a/web/content/docs/userguide/basics/project-file.md b/web/content/docs/userguide/basics/project-file.md
index dbba08b7ef7b6aa073a5941ee235316ea0707033..454cd32cf9f4ea03b44e4be576c7d53802851759 100644
--- a/web/content/docs/userguide/basics/project-file.md
+++ b/web/content/docs/userguide/basics/project-file.md
@@ -102,9 +102,29 @@ In this case you have just one patch file.
 
 ---
 
-Both methods, include and patch, can be combined.
+<div class='note'>
 
----
+### Combination of include and patch method
+
+When both methods are combined the logical order is the following:
+
+1. Apply patches
+2. Insert includes
+3. Apply patches marked with `after_includes="true"`-attribute only.
+
+Example:
+
+```xml
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="cube_1e0.prj">
+    <!-- The first line is run before the includes: -->
+    <replace sel="/*/mesh/text()">cube_1x1x1_hex_1e1.vtu</replace>
+    <!-- The following is evaluated after the includes are run: -->
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">cube_1e1</replace>
+</OpenGeoSysProjectDiff>
+```
+
+</div>
 
 ## Check project file syntax with `xmllint`