diff --git a/.gitmodules b/.gitmodules
index 9d87446c19d3de9a27f8d8d1d123db16a4ba7b00..5617fd287daa9e5cb46ee21a5d65e5ebbbffd0dc 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -46,3 +46,6 @@
 [submodule "ThirdParty/spdlog"]
 	path = ThirdParty/spdlog
 	url = https://github.com/gabime/spdlog.git
+[submodule "ThirdParty/xdmfdiff"]
+	path = ThirdParty/xdmfdiff
+	url = https://gitlab.opengeosys.org/ogs/xdmfdiff.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 900ccc90f7bcff4ee97c735b179e9fe5420e28e8..a36421bb8a009353f7113dd7305decea1e081a5a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -265,6 +265,12 @@ if(EXISTS ${PROJECT_SOURCE_DIR}/ThirdParty/vtkdiff/CMakeLists.txt AND BUILD_TEST
     install(PROGRAMS $<TARGET_FILE:vtkdiff> DESTINATION bin COMPONENT ogs_extras)
 endif()
 
+# xdmfdiff
+if(EXISTS ${PROJECT_SOURCE_DIR}/ThirdParty/xdmfdiff/CMakeLists.txt AND BUILD_TESTING)
+    add_subdirectory(ThirdParty/xdmfdiff)
+    install(PROGRAMS $<TARGET_FILE:xdmfdiff> DESTINATION bin COMPONENT ogs_extras)
+endif()
+
 include(scripts/cmake/CheckHeaderCompilation.cmake)
 
 add_subdirectory(Applications)
diff --git a/MeshLib/CMakeLists.txt b/MeshLib/CMakeLists.txt
index 572180a9cebc55df420f57dc935db2756a829f4f..da1d83bc77338536ff48aa11a59030dfe54502d6 100644
--- a/MeshLib/CMakeLists.txt
+++ b/MeshLib/CMakeLists.txt
@@ -12,6 +12,7 @@ append_source_files(SOURCES Elements)
 append_source_files(SOURCES IO)
 append_source_files(SOURCES IO/Legacy)
 append_source_files(SOURCES IO/VtkIO)
+append_source_files(SOURCES IO/XDMF)
 append_source_files(SOURCES MeshQuality)
 append_source_files(SOURCES Vtk)
 
diff --git a/MeshLib/IO/XDMF/writeXdmf.cpp b/MeshLib/IO/XDMF/writeXdmf.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa2ecdf61651ce37ade76caea78cc96f59569f75
--- /dev/null
+++ b/MeshLib/IO/XDMF/writeXdmf.cpp
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \author Tobias Meisel
+ * \date   2020-10-06
+ * \brief  Implementation of WriteXdmf3 function.
+ *
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "writeXdmf.h"
+
+#include <vtkNew.h>
+#include <vtkSmartPointer.h>
+#include <vtkXdmf3Writer.h>
+
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
+
+namespace MeshLib
+{
+namespace IO
+{
+bool writeXdmf3(const MeshLib::Mesh& mesh, std::string const& file_name)
+{
+    vtkSmartPointer<vtkXdmf3Writer> writer =
+        vtkSmartPointer<vtkXdmf3Writer>::New();
+
+    writer->SetFileName(file_name.c_str());
+    vtkNew<MeshLib::VtkMappedMeshSource> vtkSource;
+    vtkSource->SetMesh(&mesh);
+    vtkSource->Update();
+    writer->SetInputData(vtkSource->GetOutput());
+    writer->Write();
+    return true;
+}
+}  // end namespace IO
+}  // end namespace MeshLib
\ No newline at end of file
diff --git a/MeshLib/IO/XDMF/writeXdmf.h b/MeshLib/IO/XDMF/writeXdmf.h
new file mode 100644
index 0000000000000000000000000000000000000000..b723cde680734b62d18c4603253082176379eb91
--- /dev/null
+++ b/MeshLib/IO/XDMF/writeXdmf.h
@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \author Tobias Meisel
+ * \date   2020-10-06
+ * \brief  Writes MeshLib:Mesh to a Xdmf formated file
+ *
+ * \copyright
+ * Copyright (c) 2012-2020, 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 <string>
+
+namespace MeshLib
+{
+class Mesh;
+
+namespace IO
+{
+/// Writes mesh to XDMF file.
+/// \return True on success, false on error
+/// \param mesh           Mesh holds all data to be written.
+/// \param file_name      File name.
+bool writeXdmf3(MeshLib::Mesh const& mesh, std::string const& file_name);
+
+}  // end namespace IO
+}  // end namespace MeshLib
\ No newline at end of file
diff --git a/MeshLib/IO/writeMeshToFile.cpp b/MeshLib/IO/writeMeshToFile.cpp
index debc78c3bd7caee1d1aa78616e9135722b589693..256680b6a6d24f0693200aba83470483602658fc 100644
--- a/MeshLib/IO/writeMeshToFile.cpp
+++ b/MeshLib/IO/writeMeshToFile.cpp
@@ -10,14 +10,14 @@
 #include "writeMeshToFile.h"
 
 #include "BaseLib/Logging.h"
-
 #include "BaseLib/FileTools.h"
 #include "BaseLib/StringTools.h"
 
 #include "MeshLib/Mesh.h"
-
 #include "MeshLib/IO/Legacy/MeshIO.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
+#include "MeshLib/IO/XDMF/writeXdmf.h"
+
 
 namespace MeshLib
 {
@@ -44,9 +44,19 @@ int writeMeshToFile(const MeshLib::Mesh &mesh, const std::string &file_name)
         }
         return 0;
     }
-
+    if (BaseLib::hasFileExtension(".xdmf", file_name))
+    {
+        if (auto const result = writeXdmf3(mesh, file_name); !result)
+        {
+            ERR("writeMeshToFile(): Could not write mesh to '{:s}'.",
+                file_name);
+            return -1;
+        }
+        return 0;
+    }
     ERR("writeMeshToFile(): Unknown mesh file format in file {:s}.", file_name);
     return -1;
+
 }
 
 } // end namespace IO
diff --git a/MeshLib/IO/writeMeshToFile.h b/MeshLib/IO/writeMeshToFile.h
index 223a2a0d31f6f663d3ec710458d027d44859546e..d2751c7a6675f61d89623bf2e481520f29043115 100644
--- a/MeshLib/IO/writeMeshToFile.h
+++ b/MeshLib/IO/writeMeshToFile.h
@@ -18,4 +18,4 @@ namespace IO
 {
 int writeMeshToFile(const MeshLib::Mesh &mesh, const std::string &file_name);
 }
-}
+}
\ No newline at end of file
diff --git a/ProcessLib/LiquidFlow/Tests.cmake b/ProcessLib/LiquidFlow/Tests.cmake
index 1164dbcd60beade53293a971bb6c31fb0edf8517..c140ea7ae9af3c1e5ccfe974d4aeb957b17eed8a 100644
--- a/ProcessLib/LiquidFlow/Tests.cmake
+++ b/ProcessLib/LiquidFlow/Tests.cmake
@@ -430,6 +430,24 @@ AddTest(
     LF_square_1x1_tri_1.8e1_surfaceflux_ts_2_t_0.864000_expected.vtu LF_square_1x1_tri_1.8e1_surfaceflux_square_1x1_tri_1.8e1_ts_2_t_0.864000.vtu pressure pressure 1e-7 1e-13
 )
 
+AddTest(
+    NAME SimpleSynthetics_XDMF
+    PATH Parabolic/LiquidFlow/SimpleSynthetics/XDMF
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS FunctionParameterTest_XDMF.prj
+    WRAPPER time
+    TESTER xdmfdiff
+    # See https://gitlab.opengeosys.org/ogs/ogs/-/merge_requests/3184#note_85104
+    REQUIREMENTS NOT OGS_USE_MPI AND NOT COMPILER_IS_APPLE_CLANG
+    DIFF_DATA
+    square_5x5_tris_32_right_boundary_ts_0_t_0.000000.xdmf square_5x5_tris_32_right_boundary_ts_0_t_0.000000.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_right_boundary_ts_1_t_1.000000.xdmf square_5x5_tris_32_right_boundary_ts_1_t_1.000000.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_left_boundary_ts_0_t_0.000000.xdmf square_5x5_tris_32_left_boundary_ts_0_t_0.000000.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_left_boundary_ts_1_t_1.000000.xdmf square_5x5_tris_32_left_boundary_ts_1_t_1.000000.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_ts_0_t_0.000000.xdmf square_5x5_tris_32_ts_0_t_0.000000.xdmf pressure pressure 1e-7 1e-13
+    square_5x5_tris_32_ts_1_t_1.000000.xdmf square_5x5_tris_32_ts_1_t_1.000000.xdmf pressure pressure 1e-7 1e-13
+)
+
 #AddTest(
 #    NAME LiquidFlow_SimpleSynthetics_constraint_dirichlet_bc
 #    PATH Parabolic/LiquidFlow/SimpleSynthetics
diff --git a/ProcessLib/Output/CreateOutput.cpp b/ProcessLib/Output/CreateOutput.cpp
index 0d320c0b884eeb7028605f61b5f219fb774bb7b4..08a0a39f69cc81f49d9ee3b26394cd7a97c999bd 100644
--- a/ProcessLib/Output/CreateOutput.cpp
+++ b/ProcessLib/Output/CreateOutput.cpp
@@ -30,9 +30,31 @@ std::unique_ptr<Output> createOutput(
     std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
 {
     DBUG("Parse output configuration:");
-
+    OutputType const output_type = [](auto output_type)
+    {
+        try
+        {
+            const std::map<std::string, OutputType> outputType_to_enum = {
+                {"VTK", OutputType::vtk}, {"XDMF", OutputType::xdmf}};
+            auto type = outputType_to_enum.at(output_type);
+
+#ifdef USE_PETSC
+            if (type == ProcessLib::OutputType::xdmf)
+            {
+                OGS_FATAL("Parallel XDMF output is not supported (yet).");
+            }
+#endif
+
+            return type;
+        }
+        catch (std::out_of_range& e)
+        {
+            OGS_FATAL("No supported file type provided. Read `{:s}' from <output><type> \
+                in prj File. Supported: VTK, XDMF.",
+                output_type);
+        }
     //! \ogs_file_param{prj__time_loop__output__type}
-    config.checkConfigParameter("type", "VTK");
+    }(config.getConfigParameter<std::string>("type"));
 
     auto const prefix =
         //! \ogs_file_param{prj__time_loop__output__prefix}
@@ -143,7 +165,7 @@ std::unique_ptr<Output> createOutput(
         config.getConfigParameter<bool>("output_iteration_results", false);
 
     return std::make_unique<Output>(
-        output_directory, prefix, suffix, compress_output, data_mode,
+        output_directory, output_type, prefix, suffix, compress_output, data_mode,
         output_iteration_results, std::move(repeats_each_steps),
         std::move(fixed_output_times), std::move(output_data_specification),
         std::move(mesh_names_for_output), meshes);
diff --git a/ProcessLib/Output/Output.cpp b/ProcessLib/Output/Output.cpp
index aac339f160b87e4149526cb9915a4e795d34e9a9..f3663b732b44dab6a2db0f1316bc08870a864af4 100644
--- a/ProcessLib/Output/Output.cpp
+++ b/ProcessLib/Output/Output.cpp
@@ -11,12 +11,12 @@
 #include "Output.h"
 
 #include <cassert>
+#include <exception>
 #include <fstream>
 #include <vector>
 
-#include "BaseLib/Logging.h"
-
 #include "Applications/InSituLib/Adaptor.h"
+#include "BaseLib/Logging.h"
 #include "BaseLib/FileTools.h"
 #include "BaseLib/RunTime.h"
 #include "ProcessLib/Process.h"
@@ -95,9 +95,9 @@ bool Output::shallDoOutput(int timestep, double const t)
     return false;
 }
 
-Output::Output(std::string output_directory, std::string output_file_prefix,
-               std::string output_file_suffix, bool const compress_output,
-               std::string const& data_mode,
+Output::Output(std::string output_directory, OutputType output_file_type,
+               std::string output_file_prefix, std::string output_file_suffix,
+               bool const compress_output, std::string const& data_mode,
                bool const output_nonlinear_iteration_results,
                std::vector<PairRepeatEachSteps> repeats_each_steps,
                std::vector<double>&& fixed_output_times,
@@ -105,6 +105,7 @@ Output::Output(std::string output_directory, std::string output_file_prefix,
                std::vector<std::string>&& mesh_names_for_output,
                std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
     : _output_directory(std::move(output_directory)),
+      _output_file_type(output_file_type),
       _output_file_prefix(std::move(output_file_prefix)),
       _output_file_suffix(std::move(output_file_suffix)),
       _output_file_compression(compress_output),
@@ -132,12 +133,11 @@ void Output::addProcess(ProcessLib::Process const& process,
         _mesh_names_for_output.push_back(process.getMesh().getName());
     }
 
-
     for (auto const& mesh_output_name : _mesh_names_for_output)
     {
         auto const filename =
-            constructPVDName(_output_directory, _output_file_prefix,
-                             process_id, mesh_output_name);
+            constructPVDName(_output_directory, _output_file_prefix, process_id,
+                             mesh_output_name);
         _process_to_pvd_file.emplace(std::piecewise_construct,
                                      std::forward_as_tuple(&process),
                                      std::forward_as_tuple(filename));
@@ -180,16 +180,15 @@ MeshLib::IO::PVDFile* Output::findPVDFile(
 
 struct Output::OutputFile
 {
-    OutputFile(std::string const& directory, std::string const& prefix,
-               std::string const& suffix, std::string const& mesh_name,
-               int const process_id, int const timestep, double const t,
-               int const data_mode_, bool const compression_)
-        : name(BaseLib::constructFormattedFileName(prefix, mesh_name,
-                                                   process_id, timestep, t) +
-               BaseLib::constructFormattedFileName(suffix, mesh_name,
-                                                   process_id, timestep, t) +
-               ".vtu"),
+    OutputFile(std::string const& directory, OutputType const type,
+               std::string const& prefix, std::string const& suffix,
+               std::string const& mesh_name, int const process_id,
+               int const timestep, double const t, int const data_mode_,
+               bool const compression_)
+        : name(constructFilename(type, prefix, suffix, mesh_name, process_id,
+                                 timestep, t)),
           path(BaseLib::joinPaths(directory, name)),
+          type(type),
           data_mode(data_mode_),
           compression(compression_)
     {
@@ -197,6 +196,7 @@ struct Output::OutputFile
 
     std::string const name;
     std::string const path;
+    OutputType const type;
     //! Chooses vtk's data mode for output following the enumeration given
     /// in the vtkXMLWriter: {Ascii, Binary, Appended}.  See vtkXMLWriter
     /// documentation
@@ -205,19 +205,48 @@ struct Output::OutputFile
 
     //! Enables or disables zlib-compression of the output files.
     bool const compression;
+
+    static std::string constructFilename(OutputType const type, std::string prefix,
+                                         std::string suffix,
+                                         std::string mesh_name, int const process_id,
+                                         int const timestep, double const t)
+    {
+        std::map<OutputType, std::string> filetype_to_extension = {
+            {OutputType::vtk, "vtu"}, {OutputType::xdmf, "xdmf"}};
+
+        try
+        {
+            std::string extension = filetype_to_extension.at(type);
+            return BaseLib::constructFormattedFileName(
+                       prefix, mesh_name, process_id, timestep, t) +
+                   BaseLib::constructFormattedFileName(
+                       suffix, mesh_name, process_id, timestep, t) +
+                   "." + extension;
+        }
+        catch (std::out_of_range& e)
+        {
+            OGS_FATAL(
+                "No supported file type provided. Read `{:s}' from <output><type> \
+                    in prj file. Supported: VTK, XDMF.",
+                type);
+        }
+    }
 };
 
-void Output::outputBulkMesh(OutputFile const& output_file,
-                            MeshLib::IO::PVDFile* const pvd_file,
-                            MeshLib::Mesh const& mesh,
-                            double const t) const
+void Output::outputMesh(OutputFile const& output_file,
+                        MeshLib::IO::PVDFile* const pvd_file,
+                        MeshLib::Mesh const& mesh,
+                        double const t) const
 {
     DBUG("output to {:s}", output_file.path);
 
-    pvd_file->addVTUFile(output_file.name, t);
+    if (output_file.type == OutputType::vtk)
+    {
+        pvd_file->addVTUFile(output_file.name, t);
+    }
 
     makeOutput(output_file.path, mesh, output_file.compression,
-               output_file.data_mode);
+               output_file.data_mode, output_file.type);
 }
 
 void Output::doOutputAlways(Process const& process,
@@ -256,23 +285,29 @@ void Output::doOutputAlways(Process const& process,
         return;
     }
 
-    auto output_bulk_mesh = [&]() {
-        outputBulkMesh(
-            OutputFile(_output_directory, _output_file_prefix,
-                       _output_file_suffix, process.getMesh().getName(),
-                       process_id, timestep, t, _output_file_data_mode,
-                       _output_file_compression),
-            findPVDFile(process, process_id, process.getMesh().getName()),
-            process.getMesh(), t);
+    auto output_bulk_mesh = [&](MeshLib::Mesh& mesh) {
+        OutputFile const file(_output_directory, _output_file_type,
+                              _output_file_prefix, _output_file_suffix,
+                              mesh.getName(), process_id, timestep, t,
+                              _output_file_data_mode, _output_file_compression);
+
+        MeshLib::IO::PVDFile* pvd_file = nullptr;
+        if (_output_file_type == ProcessLib::OutputType::vtk)
+        {
+            pvd_file = findPVDFile(process, process_id, mesh.getName());
+        }
+        outputMesh(file, pvd_file, mesh, t);
     };
 
     for (auto const& mesh_output_name : _mesh_names_for_output)
     {
+        // process related output
         if (process.getMesh().getName() == mesh_output_name)
         {
-            output_bulk_mesh();
+            output_bulk_mesh(process.getMesh());
             continue;
         }
+        // mesh related output
         auto& mesh = *BaseLib::findElementOrError(
             begin(_meshes), end(_meshes),
             [&mesh_output_name](auto const& m) {
@@ -315,23 +350,7 @@ void Output::doOutputAlways(Process const& process,
         // output is mesh related instead of process related. This would also
         // allow for merging bulk mesh output and arbitrary mesh output.
 
-        OutputFile const output_file{_output_directory,
-                                     _output_file_prefix,
-                                     _output_file_suffix,
-                                     mesh.getName(),
-                                     process_id,
-                                     timestep,
-                                     t,
-                                     _output_file_data_mode,
-                                     _output_file_compression};
-
-        auto pvd_file = findPVDFile(process, process_id, mesh.getName());
-        pvd_file->addVTUFile(output_file.name, t);
-
-        DBUG("output to {:s}", output_file.path);
-
-        makeOutput(output_file.path, mesh, output_file.compression,
-                   output_file.data_mode);
+        output_bulk_mesh(mesh);
     }
     INFO("[time] Output of timestep {:d} took {:g} s.", timestep,
          time_output.elapsed());
@@ -410,16 +429,17 @@ void Output::doOutputNonlinearIteration(Process const& process,
     findPVDFile(process, process_id, process.getMesh().getName());
 
     std::string const output_file_name =
-        BaseLib::constructFormattedFileName(_output_file_prefix,
-                                            process.getMesh().getName(),
-                                            process_id, timestep, t) +
-        "_nliter_" + std::to_string(iteration) + ".vtu";
+        OutputFile::constructFilename(_output_file_type,_output_file_prefix,
+                                       "_ts_{:timestep}_nliter_{:time}",
+                                       process.getMesh().getName(), process_id,
+                                       timestep, iteration);
+
     std::string const output_file_path =
         BaseLib::joinPaths(_output_directory, output_file_name);
 
     DBUG("output iteration results to {:s}", output_file_path);
     makeOutput(output_file_path, process.getMesh(), _output_file_compression,
-               _output_file_data_mode);
+               _output_file_data_mode, _output_file_type);
     INFO("[time] Output took {:g} s.", time_output.elapsed());
 }
 }  // namespace ProcessLib
diff --git a/ProcessLib/Output/Output.h b/ProcessLib/Output/Output.h
index 4a14a38a9b595e853dde05b5df4def21cd0d2aff..3e7a3b4db78fbf37b128845bf08caa998f2e0df3 100644
--- a/ProcessLib/Output/Output.h
+++ b/ProcessLib/Output/Output.h
@@ -37,7 +37,8 @@ public:
     };
 
 public:
-    Output(std::string output_directory, std::string prefix, std::string suffix,
+    Output(std::string output_directory, OutputType const type,
+           std::string prefix, std::string suffix,
            bool const compress_output, std::string const& data_mode,
            bool const output_nonlinear_iteration_results,
            std::vector<PairRepeatEachSteps> repeats_each_steps,
@@ -81,13 +82,13 @@ public:
 
 private:
     struct OutputFile;
-    void outputBulkMesh(OutputFile const& output_file,
+    void outputMesh(OutputFile const& output_file,
                         MeshLib::IO::PVDFile* const pvd_file,
                         MeshLib::Mesh const& mesh,
                         double const t) const;
-
 private:
     std::string const _output_directory;
+    OutputType const _output_file_type;
     std::string const _output_file_prefix;
     std::string const _output_file_suffix;
 
diff --git a/ProcessLib/Output/ProcessOutput.cpp b/ProcessLib/Output/ProcessOutput.cpp
index ae067ab9f937d6cfee92e0218ba9a8fc2fca4a67..e1e0deeb7a1116090b74409053575238fe555454 100644
--- a/ProcessLib/Output/ProcessOutput.cpp
+++ b/ProcessLib/Output/ProcessOutput.cpp
@@ -20,6 +20,7 @@
 #include "IntegrationPointWriter.h"
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
+#include "MeshLib/IO/XDMF/writeXdmf.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
 /// Copies the ogs_version string containing the release number and the git
@@ -250,7 +251,8 @@ void addProcessDataToMesh(
 }
 
 void makeOutput(std::string const& file_name, MeshLib::Mesh const& mesh,
-                bool const compress_output, int const data_mode)
+                bool const compress_output, int const data_mode,
+                OutputType const file_type)
 {
     // Write output file
     DBUG("Writing output to '{:s}'.", file_name);
@@ -268,8 +270,20 @@ void makeOutput(std::string const& file_name, MeshLib::Mesh const& mesh,
 #endif  //_WIN32
 #endif  //__APPLE__
 
-    MeshLib::IO::VtuInterface vtu_interface(&mesh, data_mode, compress_output);
-    vtu_interface.writeToFile(file_name);
+    switch (file_type)
+    {
+        case OutputType::vtk:
+        {
+            MeshLib::IO::VtuInterface vtu_interface(&mesh, data_mode,
+                                                    compress_output);
+            vtu_interface.writeToFile(file_name);
+            break;
+        }
+        case OutputType::xdmf:
+        {
+            MeshLib::IO::writeXdmf3(mesh, file_name);
+        }
+    }
 
     // Restore floating-point exception handling.
 #ifndef _WIN32
diff --git a/ProcessLib/Output/ProcessOutput.h b/ProcessLib/Output/ProcessOutput.h
index 5e6fe96aed2c3e50a7fee8fddc5152bdb0611bf3..70b52f841ab835d969e40641762af0e80aff827c 100644
--- a/ProcessLib/Output/ProcessOutput.h
+++ b/ProcessLib/Output/ProcessOutput.h
@@ -42,11 +42,16 @@ void addProcessDataToMesh(
         integration_point_writer,
     OutputDataSpecification const& output_data_specification);
 
-//! Writes output to the given \c file_name using the VTU file format.
+//! Writes output to the given \c file_name using the specified file format.
 ///
 /// See Output::_output_file_data_mode documentation for the data_mode
 /// parameter.
+enum class OutputType : uint8_t
+{
+    vtk,
+    xdmf
+};
 void makeOutput(std::string const& file_name, MeshLib::Mesh const& mesh,
-                bool const compress_output, int const data_mode);
-
+                bool const compress_output, int const data_mode,
+                OutputType const file_type);
 }  // namespace ProcessLib
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/FunctionParameterTest_XDMF.prj b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/FunctionParameterTest_XDMF.prj
new file mode 100644
index 0000000000000000000000000000000000000000..342c735455c700deb60bc02b3db9033a5872672d
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/FunctionParameterTest_XDMF.prj
@@ -0,0 +1,186 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<!-- Almost same as FunctionParameterTest - changed XDMF output and
+boundary meshes will be written as well -->
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>square_5x5_tris_32.vtu</mesh>
+        <mesh>square_5x5_tris_32_left_boundary.vtu</mesh>
+        <mesh>square_5x5_tris_32_right_boundary.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>LiquidFlow</name>
+            <type>LIQUID_FLOW</type>
+            <integration_order>2</integration_order>
+            <darcy_gravity>
+                <axis_id>1</axis_id>
+                <g>0</g>
+            </darcy_gravity>
+            <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">
+            <phases>
+                <phase>
+                    <type>AqueousLiquid</type>
+                    <properties>
+                        <property>
+                            <name>viscosity</name>
+                            <type>Constant</type>
+                            <value> 0.0011373 </value>
+                        </property>
+                        <property>
+                            <name>density</name>
+                            <type>Constant</type>
+                            <value> 1000 </value>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+            <properties>
+                <property>
+                    <name>permeability</name>
+                    <type>Constant</type>
+                    <value>3.2439e-12</value>
+                </property>
+                <property>
+                    <name>reference_temperature</name>
+                    <type>Constant</type>
+                    <value>293.15</value>
+                </property>
+                <property>
+                    <name>porosity</name>
+                    <type>Constant</type>
+                    <value>0.17</value>
+                </property>
+                <property>
+                    <name>storage</name>
+                    <type>Constant</type>
+                    <value> 0.0 </value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <parameters>
+        <parameter>
+             <name>p_0</name>
+             <type>Function</type>
+             <expression>z</expression>
+        </parameter>
+        <parameter>
+             <name>p_polyline_left</name>
+             <type>Function</type>
+             <expression>x+1</expression>
+        </parameter>
+        <parameter>
+             <name>p_polyline_right</name>
+             <type>Function</type>
+             <expression>x-5</expression>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>pressure</name>
+            <order>1</order>
+            <components>1</components>
+            <initial_condition>p_0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <type>Dirichlet</type>
+                    <mesh>square_5x5_tris_32_left_boundary</mesh>
+                    <parameter>p_polyline_left</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <type>Dirichlet</type>
+                    <mesh>square_5x5_tris_32_right_boundary</mesh>
+                    <parameter>p_polyline_right</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <time_loop>
+        <processes>
+            <process ref="LiquidFlow">
+                <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>FixedTimeStepping</type>
+                    <t_initial> 0 </t_initial>
+                    <t_end> 1 </t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>XDMF</type>
+            <prefix>{:meshname}</prefix>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+            <meshes>
+                <mesh>square_5x5_tris_32</mesh>
+                <mesh>square_5x5_tris_32_left_boundary</mesh>
+				<mesh>square_5x5_tris_32_right_boundary</mesh>
+            </meshes>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+            <output_iteration_results>false</output_iteration_results>
+            <variables>
+                <variable> pressure </variable>
+                <variable> v </variable>
+            </variables>
+        </output>
+    </time_loop>
+    <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-15 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>1000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>lf</prefix>
+                <parameters>-lf_ksp_type cg -lf_pc_type bjacobi -lf_ksp_rtol 1e-16 -lf_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <regex>square_5x5_tris_32_ts_1_t_1.000000.vtu</regex>
+            <field>pressure</field>
+            <absolute_tolerance>1e-15</absolute_tolerance>
+            <relative_tolerance>1e-14</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32.vtu b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..af7bdb6b4b1f69d6f3aa8d7328b55cc447297416
Binary files /dev/null and b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32.vtu differ
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary.vtu b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..52b20854c8fa372bdc09ddcf39b91e323c1982ae
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary.vtu
@@ -0,0 +1,51 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="5" NumberOfCells="4">
+      <PointData Scalars="bulk_node_ids">
+        <DataArray type="UInt64" IdType="1" Name="bulk_node_ids" format="ascii" RangeMin="0" RangeMax="20">
+          5 0 10 15 20
+        </DataArray>
+      </PointData>
+      <CellData Scalars="bulk_element_ids">
+        <DataArray type="UInt64" IdType="1" Name="bulk_element_ids" format="ascii" RangeMin="0" RangeMax="24">
+          0 8 16 24
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="5">
+          0 1.25 3 0 0 3
+          0 2.5 3 0 3.75 3
+          0 5 3
+          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              5
+            </Value>
+          </InformationKey>
+          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              0
+            </Value>
+            <Value index="1">
+              5
+            </Value>
+          </InformationKey>
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="4">
+          0 1 2 0 3 2
+          4 3
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="2" RangeMax="8">
+          2 4 6 8
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="3" RangeMax="3">
+          3 3 3 3
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.h5 b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a
Binary files /dev/null and b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.h5 differ
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..4b083795314d3b95df55654400ae6c373befc54d
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_0_t_0.000000.xdmf
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="5 3" Format="XML" Precision="8">0 1.25 3 0 0 3 0 2.5 3 0 3.75 3 0 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="4" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="16" Format="XML" Precision="8">2 2 0 1 2 2 2 0 2 2 3 2 2 2 4 3</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 52 50 56 45 103 57 98 49 100 97 48 98 56 57 46 100 105 114 116 121 46 50 48 50 48 49 48 49 50 48 55 52 56 52 52</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="5" Format="XML" Precision="8">3 3 3 3 3</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.h5 b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a
Binary files /dev/null and b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.h5 differ
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..16cb9288e16d27dea51a00203ea718a4add28889
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_left_boundary_ts_1_t_1.000000.xdmf
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="5 3" Format="XML" Precision="8">0 1.25 3 0 0 3 0 2.5 3 0 3.75 3 0 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="4" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="16" Format="XML" Precision="8">2 2 0 1 2 2 2 0 2 2 3 2 2 2 4 3</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 52 50 56 45 103 57 98 49 100 97 48 98 56 57 46 100 105 114 116 121 46 50 48 50 48 49 48 49 50 48 55 52 56 52 52</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="5" Format="XML" Precision="8">1 1 1 1 1</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary.vtu b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8df915d4a239106d1144e779cc796a10dd270e87
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary.vtu
@@ -0,0 +1,51 @@
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="5" NumberOfCells="4">
+      <PointData Scalars="bulk_node_ids">
+        <DataArray type="UInt64" IdType="1" Name="bulk_node_ids" format="ascii" RangeMin="4" RangeMax="24">
+          4 9 14 19 24
+        </DataArray>
+      </PointData>
+      <CellData Scalars="bulk_element_ids">
+        <DataArray type="UInt64" IdType="1" Name="bulk_element_ids" format="ascii" RangeMin="7" RangeMax="31">
+          7 15 23 31
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="5" RangeMax="7.0710678118654755">
+          5 0 3 5 1.25 3
+          5 2.5 3 5 3.75 3
+          5 5 3
+          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              5
+            </Value>
+            <Value index="1">
+              7.0710678119
+            </Value>
+          </InformationKey>
+          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
+            <Value index="0">
+              5
+            </Value>
+            <Value index="1">
+              7.0710678119
+            </Value>
+          </InformationKey>
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="4">
+          0 1 1 2 2 3
+          3 4
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="2" RangeMax="8">
+          2 4 6 8
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="3" RangeMax="3">
+          3 3 3 3
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.h5 b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a
Binary files /dev/null and b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.h5 differ
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..19e83b7ee362229dcdc347096cf372490bbc6eaf
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_0_t_0.000000.xdmf
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="5 3" Format="XML" Precision="8">5 0 3 5 1.25 3 5 2.5 3 5 3.75 3 5 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="4" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="16" Format="XML" Precision="8">2 2 0 1 2 2 1 2 2 2 2 3 2 2 3 4</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 52 50 56 45 103 57 98 49 100 97 48 98 56 57 46 100 105 114 116 121 46 50 48 50 48 49 48 49 50 48 55 52 56 52 52</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="5" Format="XML" Precision="8">3 3 3 3 3</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.h5 b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a
Binary files /dev/null and b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.h5 differ
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..76fd582bb3203c3963ddc6972b8a0b4efc2864a6
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_right_boundary_ts_1_t_1.000000.xdmf
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="5 3" Format="XML" Precision="8">5 0 3 5 1.25 3 5 2.5 3 5 3.75 3 5 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="4" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="16" Format="XML" Precision="8">2 2 0 1 2 2 1 2 2 2 2 3 2 2 3 4</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 52 50 56 45 103 57 98 49 100 97 48 98 56 57 46 100 105 114 116 121 46 50 48 50 48 49 48 49 50 48 55 52 56 52 52</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="5" Format="XML" Precision="8">0 0 0 0 0</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_0_t_0.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_0_t_0.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..489226749d6631690c1168a744567424c84ac5b6
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_0_t_0.000000.xdmf
@@ -0,0 +1,28 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="25 3" Format="XML" Precision="8">0 0 3 1.25 0 3 2.5 0 3 3.75 0 3 5 0 3 0 1.25 3 1.25 1.25 3 2.5 1.25 3 3.75 1.25 3 5 1.25 3 0 2.5 3 1.25 2.5 3 2.5 2.5 3 3.75 2.5 3 5 2.5 3 0 3.75 3 1.25 3.75 3 2.5 3.75 3 3.75 3.75 3 5 3.75 3 0 5 3 1.25 5 3 2.5 5 3 3.75 5 3 5 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="32" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="128" Format="HDF" Precision="8">square_5x5_tris_32_ts_0_t_0.000000.h5:Data7</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 53 56 54 45 103 56 102 56 53 56 53 54 49 97 46 100 105 114 116 121 46 50 48 50 48 49 48 50 54 49 51 50 52 51 53</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="HydraulicFlow" Type="None">
+        <DataItem DataType="Float" Dimensions="25" Format="XML" Precision="8">0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="MaterialIDs" Type="None">
+        <DataItem DataType="Int" Dimensions="25" Format="XML" Precision="4">0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="25" Format="XML" Precision="8">3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="v" Type="None">
+        <DataItem DataType="Float" Dimensions="25 2" Format="XML" Precision="8">0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_1_t_1.000000.xdmf b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_1_t_1.000000.xdmf
new file mode 100644
index 0000000000000000000000000000000000000000..e9a508ac22ff4ffbca83868dde7b34bf96e2a8ca
--- /dev/null
+++ b/Tests/Data/Parabolic/LiquidFlow/SimpleSynthetics/XDMF/square_5x5_tris_32_ts_1_t_1.000000.xdmf
@@ -0,0 +1,28 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
+  <Domain>
+    <Grid Name="Grid">
+      <Geometry Origin="" Type="XYZ">
+        <DataItem DataType="Float" Dimensions="25 3" Format="XML" Precision="8">0 0 3 1.25 0 3 2.5 0 3 3.75 0 3 5 0 3 0 1.25 3 1.25 1.25 3 2.5 1.25 3 3.75 1.25 3 5 1.25 3 0 2.5 3 1.25 2.5 3 2.5 2.5 3 3.75 2.5 3 5 2.5 3 0 3.75 3 1.25 3.75 3 2.5 3.75 3 3.75 3.75 3 5 3.75 3 0 5 3 1.25 5 3 2.5 5 3 3.75 5 3 5 5 3</DataItem>
+      </Geometry>
+      <Topology Dimensions="32" Type="Mixed">
+        <DataItem DataType="Int" Dimensions="128" Format="HDF" Precision="8">square_5x5_tris_32_ts_1_t_1.000000.h5:Data7</DataItem>
+      </Topology>
+      <Attribute Center="Grid" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="OGS_VERSION" Type="None">
+        <DataItem DataType="Char" Dimensions="41" Format="XML" Precision="1">54 46 51 46 50 45 53 56 54 45 103 56 102 56 53 56 53 54 49 97 46 100 105 114 116 121 46 50 48 50 48 49 48 50 54 49 51 50 52 51 53</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="HydraulicFlow" Type="None">
+        <DataItem DataType="Float" Dimensions="25" Format="XML" Precision="8">-2.8522817198628322e-09 2.8522817198628317e-09 0 4.2784225797942485e-09 -4.2784225797942485e-09 -5.7045634397256652e-09 5.7045634397256635e-09 0 8.5568451595884969e-09 -8.5568451595884969e-09 -5.7045634397256652e-09 5.7045634397256635e-09 0 8.5568451595884969e-09 -8.5568451595884969e-09 -5.7045634397256652e-09 5.7045634397256635e-09 0 8.5568451595884969e-09 -8.5568451595884969e-09 -2.8522817198628322e-09 2.8522817198628317e-09 0 4.2784225797942485e-09 -4.2784225797942485e-09</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="MaterialIDs" Type="None">
+        <DataItem DataType="Int" Dimensions="25" Format="XML" Precision="4">0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="pressure" Type="None">
+        <DataItem DataType="Float" Dimensions="25" Format="XML" Precision="8">1 0.75000000000000011 0.50000000000000067 0.25000000000000044 0 1 0.75000000000000033 0.50000000000000067 0.25000000000000044 0 1 0.75 0.50000000000000067 0.2500000000000005 0 1 0.75000000000000033 0.50000000000000078 0.2500000000000005 0 1 0.75000000000000011 0.50000000000000078 0.2500000000000005 0</DataItem>
+      </Attribute>
+      <Attribute Center="Node" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="v" Type="None">
+        <DataItem DataType="Float" Dimensions="25 2" Format="XML" Precision="8">5.7045634397256606e-10 -1.7756826847496619e-25 5.7045634397256585e-10 -2.3675769129995494e-25 5.7045634397256668e-10 0 5.7045634397256761e-10 0 5.7045634397256771e-10 0 5.7045634397256616e-10 2.3675769129995494e-25 5.7045634397256596e-10 1.1837884564997749e-25 5.7045634397256637e-10 -8.8784134237483108e-26 5.7045634397256751e-10 -5.9189422824988735e-26 5.7045634397256771e-10 0 5.7045634397256637e-10 -2.3675769129995494e-25 5.7045634397256585e-10 -2.9594711412494396e-26 5.7045634397256616e-10 5.9189422824988746e-26 5.7045634397256751e-10 -5.9189422824988735e-26 5.7045634397256792e-10 -5.9189422824988746e-26 5.7045634397256596e-10 1.1837884564997747e-25 5.7045634397256596e-10 -1.1837884564997742e-25 5.7045634397256637e-10 -1.7756826847496622e-25 5.7045634397256761e-10 -2.9594711412494373e-26 5.7045634397256792e-10 0 5.7045634397256637e-10 0 5.7045634397256585e-10 2.3675769129995494e-25 5.7045634397256596e-10 1.1837884564997749e-25 5.7045634397256751e-10 0 5.7045634397256792e-10 0</DataItem>
+      </Attribute>
+    </Grid>
+  </Domain>
+</Xdmf>
diff --git a/ThirdParty/container-maker b/ThirdParty/container-maker
index ada9862a5db500ed31b2c80c3d660ccdb6c016c3..a24d8405cdbb64e15d78a3f21f64c8fc84b82f76 160000
--- a/ThirdParty/container-maker
+++ b/ThirdParty/container-maker
@@ -1 +1 @@
-Subproject commit ada9862a5db500ed31b2c80c3d660ccdb6c016c3
+Subproject commit a24d8405cdbb64e15d78a3f21f64c8fc84b82f76
diff --git a/ThirdParty/xdmfdiff b/ThirdParty/xdmfdiff
new file mode 160000
index 0000000000000000000000000000000000000000..428760e717acef51fdedcc2bf3e56474c1b971ef
--- /dev/null
+++ b/ThirdParty/xdmfdiff
@@ -0,0 +1 @@
+Subproject commit 428760e717acef51fdedcc2bf3e56474c1b971ef
diff --git a/scripts/cmake/ConanSetup.cmake b/scripts/cmake/ConanSetup.cmake
index 260d62a557a75fa56861ff41f79b225fb7014ce2..534bc04ac736c8f27c1c45317e885bc0d5c4d689 100644
--- a/scripts/cmake/ConanSetup.cmake
+++ b/scripts/cmake/ConanSetup.cmake
@@ -41,6 +41,7 @@ set(CONAN_OPTIONS
     boost:header_only=True
     vtk:minimal=True
     vtk:ioxml=True
+    vtk:ioxdmf3=True
     CACHE INTERNAL ""
 )
 
diff --git a/scripts/cmake/Find.cmake b/scripts/cmake/Find.cmake
index 14babf80932df2b555ed6a1b11e0cbbe0a5a9903..e85d3a0ae5117b17ee09e335a5f1763a24fc3a40 100644
--- a/scripts/cmake/Find.cmake
+++ b/scripts/cmake/Find.cmake
@@ -62,7 +62,7 @@ find_program(PARSL parsl-visualize HINTS ${LOCAL_VIRTUALENV_BIN_DIRS})
 ######################
 find_package(Boost ${ogs.minimum_version.boost} REQUIRED)
 
-set(VTK_COMPONENTS vtkIOXML)
+set(VTK_COMPONENTS vtkIOXML vtkIOXdmf3)
 if(OGS_BUILD_GUI)
     set(VTK_COMPONENTS ${VTK_COMPONENTS}
         vtkIOLegacy vtkIOExport vtkImagingCore
@@ -80,6 +80,12 @@ else()
     find_package(VTK ${ogs.minimum_version.vtk} REQUIRED COMPONENTS ${VTK_COMPONENTS})
     include(${VTK_USE_FILE})
 endif()
+include_directories(SYSTEM
+    # Xdmf:
+    ${VTK_INSTALL_PREFIX}/include
+    # libxml:
+    ${VTK_INSTALL_PREFIX}/include/vtk-${VTK_VERSION_MAJOR}.${VTK_VERSION_MINOR}/vtklibxml2
+)
 
 find_package(Eigen3 ${ogs.minimum_version.eigen} REQUIRED)
 include_directories(SYSTEM ${EIGEN3_INCLUDE_DIR})
diff --git a/scripts/cmake/SubmoduleSetup.cmake b/scripts/cmake/SubmoduleSetup.cmake
index 11990bac6b6c45f63e07c23fa5def2a559263be5..7db7db3685b2b303b3566faf131904af887ac255 100644
--- a/scripts/cmake/SubmoduleSetup.cmake
+++ b/scripts/cmake/SubmoduleSetup.cmake
@@ -21,7 +21,7 @@ set(REQUIRED_SUBMODULES
     ${OGS_ADDITIONAL_SUBMODULES_TO_CHECKOUT}
 )
 if(BUILD_TESTING)
-    list(APPEND REQUIRED_SUBMODULES ThirdParty/vtkdiff)
+    list(APPEND REQUIRED_SUBMODULES ThirdParty/vtkdiff ThirdParty/xdmfdiff)
 endif()
 if(OGS_BUILD_UTILS)
     # Required by the partmesh tool, which is build with utils only.
diff --git a/scripts/cmake/test/AddTest.cmake b/scripts/cmake/test/AddTest.cmake
index e44273373efe8ce411416ea4e44db150377d8864..08dfceda927b57193733d9a02530994796ae57ca 100644
--- a/scripts/cmake/test/AddTest.cmake
+++ b/scripts/cmake/test/AddTest.cmake
@@ -136,6 +136,9 @@ function (AddTest)
     if(AddTest_TESTER STREQUAL "vtkdiff" AND NOT TARGET vtkdiff)
         return()
     endif()
+    if(AddTest_TESTER STREQUAL "xdmfdiff" AND NOT TARGET xdmfdiff)
+        return()
+    endif()
     if(AddTest_TESTER STREQUAL "memcheck" AND NOT GREP_TOOL_PATH)
         return()
     endif()
@@ -147,7 +150,7 @@ function (AddTest)
         endif()
     endif()
 
-    if((AddTest_TESTER STREQUAL "diff" OR AddTest_TESTER STREQUAL "vtkdiff") AND NOT AddTest_DIFF_DATA)
+    if((AddTest_TESTER STREQUAL "diff" OR AddTest_TESTER STREQUAL "vtkdiff" OR AddTest_TESTER STREQUAL "xdmfdiff") AND NOT AddTest_DIFF_DATA)
         message(FATAL_ERROR "AddTest(): ${AddTest_NAME} - no DIFF_DATA given!")
     endif()
 
@@ -156,6 +159,8 @@ function (AddTest)
         set(TESTER_ARGS "-sbB")
     elseif(AddTest_TESTER STREQUAL "vtkdiff")
         set(SELECTED_DIFF_TOOL_PATH $<TARGET_FILE:vtkdiff>)
+    elseif(AddTest_TESTER STREQUAL "xdmfdiff")
+        set(SELECTED_DIFF_TOOL_PATH $<TARGET_FILE:xdmfdiff>)
     endif()
 
     if(AddTest_TESTER STREQUAL "diff")
@@ -165,7 +170,7 @@ function (AddTest)
                 ${TESTER_ARGS} ${AddTest_SOURCE_PATH}/${FILE_EXPECTED} \
                 ${AddTest_BINARY_PATH}/${FILE}")
         endforeach()
-    elseif(AddTest_TESTER STREQUAL "vtkdiff")
+    elseif(AddTest_TESTER STREQUAL "vtkdiff" OR AddTest_TESTER STREQUAL "xdmfdiff")
         list(LENGTH AddTest_DIFF_DATA DiffDataLength)
         math(EXPR DiffDataLengthMod4 "${DiffDataLength} % 4")
         math(EXPR DiffDataLengthMod6 "${DiffDataLength} % 6")