diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 07a48a655155ae4316995fe7a4e82cdb42423412..e5176f05673f6c9be2ba6e8ded82e8192ed08473 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -23,16 +23,13 @@
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MeshLib/Mesh.h"
 
-#include "NumLib/ODESolver/TimeDiscretizationBuilder.h"
 #include "NumLib/ODESolver/ConvergenceCriterion.h"
 
 // FileIO
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "MeshLib/IO/readMeshFromFile.h"
 
-#include "BaseLib/ConfigTree.h"
-
-#include "UncoupledProcessesTimeLoop.h"
+#include "ProcessLib/UncoupledProcessesTimeLoop.h"
 
 #include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
@@ -85,17 +82,14 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
     //! \ogs_file_param{prj__processes}
     parseProcesses(project_config.getConfigSubtree("processes"));
 
-    //! \ogs_file_param{prj__output}
-    parseOutput(project_config.getConfigSubtree("output"), output_directory);
-
-    //! \ogs_file_param{prj__time_stepping}
-    parseTimeStepping(project_config.getConfigSubtree("time_stepping"));
-
     //! \ogs_file_param{prj__linear_solvers}
     parseLinearSolvers(project_config.getConfigSubtree("linear_solvers"));
 
     //! \ogs_file_param{prj__nonlinear_solvers}
     parseNonlinearSolvers(project_config.getConfigSubtree("nonlinear_solvers"));
+
+    //! \ogs_file_param{prj__time_loop}
+    parseTimeLoop(project_config.getConfigSubtree("time_loop"), output_directory);
 }
 
 ProjectData::~ProjectData()
@@ -152,57 +146,6 @@ bool ProjectData::removeMesh(const std::string& name)
     return mesh_found;
 }
 
-void ProjectData::buildProcesses()
-{
-    for (auto const& pc : _process_configs)
-    {
-        //! \ogs_file_param{process__type}
-        auto const type = pc.peekConfigParameter<std::string>("type");
-
-        auto const nl_slv_name =
-            //! \ogs_file_param{process__nonlinear_solver}
-            pc.getConfigParameter<std::string>("nonlinear_solver");
-        auto& nl_slv = BaseLib::getOrError(
-            _nonlinear_solvers, nl_slv_name,
-            "A nonlinear solver with the given name has not been defined.");
-
-        auto time_disc = NumLib::createTimeDiscretization(
-            //! \ogs_file_param{process__time_discretization}
-            pc.getConfigSubtree("time_discretization"));
-
-        auto conv_crit = NumLib::createConvergenceCriterion(
-            //! \ogs_file_param{process__convergence_criterion}
-            pc.getConfigSubtree("convergence_criterion"));
-
-        if (type == "GROUNDWATER_FLOW")
-        {
-            // The existence check of the in the configuration referenced
-            // process variables is checked in the physical process.
-            // TODO at the moment we have only one mesh, later there can be
-            // several meshes. Then we have to assign the referenced mesh
-            // here.
-            _processes.emplace_back(
-                ProcessLib::GroundwaterFlow::createGroundwaterFlowProcess(
-                    *_mesh_vec[0], *nl_slv, std::move(time_disc),
-                    std::move(conv_crit), _process_variables, _parameters, pc));
-        }
-        else if (type == "TES")
-        {
-            _processes.emplace_back(ProcessLib::TES::createTESProcess(
-                *_mesh_vec[0], *nl_slv, std::move(time_disc),
-                std::move(conv_crit), _process_variables, _parameters, pc));
-        }
-        else
-        {
-            OGS_FATAL("Unknown process type: %s", type.c_str());
-        }
-    }
-
-    // process configs are not needed anymore, so clear the storage
-    // in order to trigger config tree checks
-    _process_configs.clear();
-}
-
 bool ProjectData::meshExists(const std::string& name) const
 {
     return findMeshByName(name) != _mesh_vec.end();
@@ -297,32 +240,48 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)
     //! \ogs_file_param{prj__processes__process}
     for (auto process_config : processes_config.getConfigSubtreeList("process"))
     {
-        // process type must be specified.
         //! \ogs_file_param{process__type}
-        process_config.peekConfigParameter<std::string>("type");
-        process_config.ignoreConfigParameter("name");
-        _process_configs.push_back(std::move(process_config));
-    }
-}
+        auto const type = process_config.peekConfigParameter<std::string>("type");
 
-void ProjectData::parseOutput(BaseLib::ConfigTree const& output_config,
-                              std::string const& output_directory)
-{
-    //! \ogs_file_param_special{prj__output__VTK}
-    //! \ogs_file_param{prj__output__type}
-    output_config.checkConfigParameter("type", "VTK");
-    DBUG("Parse output configuration:");
+        //! \ogs_file_param{process__type}
+        auto const name = process_config.getConfigParameter<std::string>("name");
+
+        std::unique_ptr<ProcessLib::Process> process;
+
+        if (type == "GROUNDWATER_FLOW")
+        {
+            // The existence check of the in the configuration referenced
+            // process variables is checked in the physical process.
+            // TODO at the moment we have only one mesh, later there can be
+            // several meshes. Then we have to assign the referenced mesh
+            // here.
+            process = ProcessLib::GroundwaterFlow::createGroundwaterFlowProcess(
+                *_mesh_vec[0], _process_variables, _parameters, process_config);
+        }
+        else if (type == "TES")
+        {
+            process = ProcessLib::TES::createTESProcess(
+                *_mesh_vec[0], _process_variables, _parameters, process_config);
+        }
+        else
+        {
+            OGS_FATAL("Unknown process type: %s", type.c_str());
+        }
 
-    _output = ProcessLib::Output::newInstance(output_config, output_directory);
+        BaseLib::insertIfKeyUniqueElseError(_processes,
+                                            name,
+                                            std::move(process),
+                                            "The process name is not unique");
+    }
 }
 
-void ProjectData::parseTimeStepping(
-    BaseLib::ConfigTree const& timestepping_config)
+void ProjectData::parseTimeLoop(BaseLib::ConfigTree const& config,
+                                std::string const& output_directory)
 {
     DBUG("Reading time loop configuration.");
 
-    _time_loop =
-        ApplicationsLib::createUncoupledProcessesTimeLoop(timestepping_config);
+    _time_loop = ProcessLib::createUncoupledProcessesTimeLoop(
+        config, output_directory, _processes, _nonlinear_solvers);
 
     if (!_time_loop)
     {
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index ed8006bfd1767858f4b71af1b6b1d8f1cceb373b..8bed77758b1ca4db94bf8ad10213917e3193bb5d 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -22,13 +22,10 @@
 
 #include "GeoLib/GEOObjects.h"
 
-#include "ProcessLib/Output.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/ProcessVariable.h"
 
-#include "NumLib/ODESolver/Types.h"
-
 namespace MathLib
 {
 class PiecewiseLinearInterpolation;
@@ -43,7 +40,7 @@ namespace NumLib
 class NonlinearSolverBase;
 }
 
-namespace ApplicationsLib
+namespace ProcessLib
 {
 class UncoupledProcessesTimeLoop;
 }
@@ -57,7 +54,7 @@ class ProjectData final
 {
 public:
     /// The time loop type used to solve this project's processes.
-    using TimeLoop = ApplicationsLib::UncoupledProcessesTimeLoop;
+    using TimeLoop = ProcessLib::UncoupledProcessesTimeLoop;
 
     /// The empty constructor used in the gui, for example, when the project's
     /// configuration is not loaded yet.
@@ -105,35 +102,15 @@ public:
     // Process interface
     //
 
-    /// Builds processes.
-    void buildProcesses();
-
-    /// Iterator access for processes.
     /// Provides read access to the process container.
-    std::vector<std::unique_ptr<ProcessLib::Process>>::const_iterator
-    processesBegin() const
-    {
-        return _processes.begin();
-    }
-    std::vector<std::unique_ptr<ProcessLib::Process>>::iterator processesBegin()
+    std::map<std::string, std::unique_ptr<ProcessLib::Process>> const&
+    getProcesses() const
     {
-        return _processes.begin();
+        return _processes;
     }
 
-    /// Iterator access for processes as in processesBegin().
-    std::vector<std::unique_ptr<ProcessLib::Process>>::const_iterator
-    processesEnd() const
-    {
-        return _processes.end();
-    }
-    std::vector<std::unique_ptr<ProcessLib::Process>>::iterator processesEnd()
-    {
-        return _processes.end();
-    }
-
-    ProcessLib::Output const& getOutputControl() const { return *_output; }
-    ProcessLib::Output& getOutputControl() { return *_output; }
     TimeLoop& getTimeLoop() { return *_time_loop; }
+
 private:
     /// Checks if a mesh with the same name exists and provides a unique name in
     /// case of already existing mesh. Returns true if the mesh name is unique.
@@ -164,12 +141,9 @@ private:
     /// constructor.
     void parseProcesses(BaseLib::ConfigTree const& process_config);
 
-    /// Parses the output configuration.
-    /// Parses the file tag and sets output file prefix.
-    void parseOutput(BaseLib::ConfigTree const& output_config,
-                     std::string const& output_directory);
-
-    void parseTimeStepping(BaseLib::ConfigTree const& timestepping_config);
+    /// Parses the time loop configuration.
+    void parseTimeLoop(BaseLib::ConfigTree const& config,
+                       const std::string& output_directory);
 
     void parseLinearSolvers(BaseLib::ConfigTree const& config);
 
@@ -177,29 +151,21 @@ private:
 
     void parseCurves(boost::optional<BaseLib::ConfigTree> const& config);
 
-private:
     GeoLib::GEOObjects* _geoObjects = new GeoLib::GEOObjects();
     std::vector<MeshLib::Mesh*> _mesh_vec;
-    std::vector<std::unique_ptr<ProcessLib::Process>> _processes;
+    std::map<std::string, std::unique_ptr<ProcessLib::Process>> _processes;
     std::vector<ProcessLib::ProcessVariable> _process_variables;
 
-    /// Buffer for each process' config used in the process building function.
-    std::vector<BaseLib::ConfigTree> _process_configs;
-
     /// Buffer for each parameter config passed to the process.
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> _parameters;
 
-    std::unique_ptr<ProcessLib::Output> _output;
-
     /// The time loop used to solve this project's processes.
     std::unique_ptr<TimeLoop> _time_loop;
 
-    std::map<std::string,
-             std::unique_ptr<GlobalLinearSolver>>
-        _linear_solvers;
+    std::map<std::string, std::unique_ptr<GlobalLinearSolver>> _linear_solvers;
 
-    using NonlinearSolver = NumLib::NonlinearSolverBase;
-    std::map<std::string, std::unique_ptr<NonlinearSolver>> _nonlinear_solvers;
+    std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>>
+        _nonlinear_solvers;
     std::map<std::string,
              std::unique_ptr<MathLib::PiecewiseLinearInterpolation>>
         _curves;
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
deleted file mode 100644
index cdfdecec8e2740576e9c07d3caf3b9a306669d63..0000000000000000000000000000000000000000
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
+++ /dev/null
@@ -1,257 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, 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 "UncoupledProcessesTimeLoop.h"
-#include "BaseLib/RunTime.h"
-
-namespace ApplicationsLib
-{
-std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
-    BaseLib::ConfigTree const& conf)
-{
-    //! \ogs_file_param{prj__time_stepping__type}
-    auto const type = conf.peekConfigParameter<std::string>("type");
-
-    std::unique_ptr<NumLib::ITimeStepAlgorithm> timestepper;
-
-    if (type == "SingleStep")
-    {
-        //! \ogs_file_param_special{prj__time_stepping__SingleStep}
-        conf.ignoreConfigParameter("type");
-        timestepper.reset(new NumLib::FixedTimeStepping(0.0, 1.0, 1.0));
-    }
-    else if (type == "FixedTimeStepping")
-    {
-        timestepper = NumLib::FixedTimeStepping::newInstance(conf);
-    }
-    else
-    {
-        OGS_FATAL("Unknown timestepper type: `%s'.", type.c_str());
-    }
-
-    using TimeLoop = UncoupledProcessesTimeLoop;
-    return std::unique_ptr<TimeLoop>{new TimeLoop{std::move(timestepper)}};
-}
-
-std::vector<typename UncoupledProcessesTimeLoop::SingleProcessData>
-UncoupledProcessesTimeLoop::initInternalData(ProjectData& project)
-{
-    auto const num_processes =
-        std::distance(project.processesBegin(), project.processesEnd());
-
-    std::vector<SingleProcessData> per_process_data;
-    per_process_data.reserve(num_processes);
-
-    // create a time discretized ODE system for each process
-    for (auto p = project.processesBegin(); p != project.processesEnd(); ++p)
-    {
-        auto& pcs = **p;
-        auto& nonlinear_solver = pcs.getNonlinearSolver();
-        auto& time_disc = pcs.getTimeDiscretization();
-
-        per_process_data.emplace_back(
-            makeSingleProcessData(nonlinear_solver, **p, time_disc));
-    }
-
-    return per_process_data;
-}
-
-void UncoupledProcessesTimeLoop::setInitialConditions(
-    ProjectData& project,
-    double const t0,
-    std::vector<SingleProcessData>& per_process_data)
-{
-    auto const num_processes =
-        std::distance(project.processesBegin(), project.processesEnd());
-
-    _process_solutions.reserve(num_processes);
-
-    unsigned pcs_idx = 0;
-    for (auto p = project.processesBegin(); p != project.processesEnd();
-         ++p, ++pcs_idx)
-    {
-        auto& pcs = **p;
-        auto& time_disc = pcs.getTimeDiscretization();
-
-        auto& ppd = per_process_data[pcs_idx];
-        auto& ode_sys = *ppd.tdisc_ode_sys;
-        auto const nl_tag = ppd.nonlinear_solver_tag;
-
-        // append a solution vector of suitable size
-        _process_solutions.emplace_back(
-            &NumLib::GlobalVectorProvider::provider.getVector(
-                ode_sys.getMatrixSpecifications()));
-
-        auto& x0 = *_process_solutions[pcs_idx];
-        pcs.setInitialConditions(t0, x0);
-        MathLib::LinAlg::finalizeAssembly(x0);
-
-        time_disc.setInitialState(t0, x0);  // push IC
-
-        if (time_disc.needsPreload())
-        {
-            auto& nonlinear_solver = ppd.nonlinear_solver;
-            auto& mat_strg = ppd.mat_strg;
-            auto& conv_crit = pcs.getConvergenceCriterion();
-
-            setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
-            nonlinear_solver.assemble(x0);
-            time_disc.pushState(
-                t0, x0, mat_strg);  // TODO: that might do duplicate work
-        }
-    }
-}
-
-
-bool
-UncoupledProcessesTimeLoop::
-solveOneTimeStepOneProcess(
-        GlobalVector& x, std::size_t const timestep, double const t, double const delta_t,
-        SingleProcessData& process_data,
-        UncoupledProcessesTimeLoop::Process& process,
-        ProcessLib::Output const& output_control)
-{
-    auto& time_disc = process.getTimeDiscretization();
-    auto& conv_crit = process.getConvergenceCriterion();
-    auto& ode_sys = *process_data.tdisc_ode_sys;
-    auto& nonlinear_solver = process_data.nonlinear_solver;
-    auto const nl_tag = process_data.nonlinear_solver_tag;
-
-    setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
-
-    // Note: Order matters!
-    // First advance to the next timestep, then set known solutions at that
-    // time, afterwards pass the right solution vector and time to the
-    // preTimestep() hook.
-
-    time_disc.nextTimestep(t, delta_t);
-
-    applyKnownSolutions(ode_sys, nl_tag, x);
-
-    process.preTimestep(x, t, delta_t);
-
-    auto const post_iteration_callback = [&](unsigned iteration,
-                                             GlobalVector const& x) {
-        output_control.doOutputNonlinearIteration(process, timestep, t, x, iteration);
-    };
-
-    bool nonlinear_solver_succeeded =
-        nonlinear_solver.solve(x, post_iteration_callback);
-
-    auto& mat_strg = process_data.mat_strg;
-    time_disc.pushState(t, x, mat_strg);
-
-    process.postTimestep(x);
-
-    return nonlinear_solver_succeeded;
-}
-
-bool UncoupledProcessesTimeLoop::loop(ProjectData& project)
-{
-    auto per_process_data = initInternalData(project);
-
-    auto& out_ctrl = project.getOutputControl();
-    out_ctrl.initialize(project.processesBegin(), project.processesEnd());
-
-    auto const t0 = _timestepper->getTimeStep().current();  // time of the IC
-
-    // init solution storage
-    setInitialConditions(project, t0, per_process_data);
-
-    // output initial conditions
-    {
-        unsigned pcs_idx = 0;
-        for (auto p = project.processesBegin(); p != project.processesEnd();
-             ++p, ++pcs_idx)
-        {
-            auto const& x0 = *_process_solutions[pcs_idx];
-            out_ctrl.doOutput(**p, 0, t0, x0);
-        }
-    }
-
-    double t = t0;
-    std::size_t timestep = 1;  // the first timestep really is number one
-    bool nonlinear_solver_succeeded = true;
-
-    while (_timestepper->next())
-    {
-        BaseLib::RunTime time_timestep;
-        time_timestep.start();
-
-        auto const ts = _timestepper->getTimeStep();
-        auto const delta_t = ts.dt();
-        t = ts.current();
-        timestep = ts.steps();
-
-        INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
-             timestep, t, delta_t);
-
-        // TODO use process name
-        unsigned pcs_idx = 0;
-        for (auto p = project.processesBegin(); p != project.processesEnd();
-             ++p, ++pcs_idx)
-        {
-            BaseLib::RunTime time_timestep_process;
-            time_timestep_process.start();
-
-            auto& x = *_process_solutions[pcs_idx];
-
-            nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
-                x, timestep, t, delta_t, per_process_data[pcs_idx], **p,
-                out_ctrl);
-
-            INFO("[time] Solving process #%u took %g s in timestep #%u.",
-                 timestep, time_timestep.elapsed(), pcs_idx);
-
-            if (!nonlinear_solver_succeeded)
-            {
-                ERR("The nonlinear solver failed in timestep #%u at t = %g s"
-                    " for process #%u.",
-                    timestep, t, pcs_idx);
-
-                // save unsuccessful solution
-                out_ctrl.doOutputAlways(**p, timestep, t, x);
-
-                break;
-            }
-            else
-            {
-                out_ctrl.doOutput(**p, timestep, t, x);
-            }
-        }
-
-        INFO("[time] Timestep #%u took %g s.", timestep,
-             time_timestep.elapsed());
-
-        if (!nonlinear_solver_succeeded)
-            break;
-    }
-
-    // output last timestep
-    if (nonlinear_solver_succeeded)
-    {
-        unsigned pcs_idx = 0;
-        for (auto p = project.processesBegin(); p != project.processesEnd();
-             ++p, ++pcs_idx)
-        {
-            auto const& x = *_process_solutions[pcs_idx];
-            out_ctrl.doOutputLastTimestep(**p, timestep, t, x);
-        }
-    }
-
-    return nonlinear_solver_succeeded;
-}
-
-UncoupledProcessesTimeLoop::~UncoupledProcessesTimeLoop()
-{
-    for (auto* x : _process_solutions)
-        NumLib::GlobalVectorProvider::provider.releaseVector(*x);
-}
-
-}  // namespace ApplicationsLib
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
deleted file mode 100644
index 7cdce87b19e3bdf3174e71552567ca30f67d6e1e..0000000000000000000000000000000000000000
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
+++ /dev/null
@@ -1,239 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef APPLICATIONSLIB_UNCOUPLED_PROCESSES_TIMELOOP
-#define APPLICATIONSLIB_UNCOUPLED_PROCESSES_TIMELOOP
-
-#include <memory>
-
-#include <logog/include/logog.hpp>
-
-#include "BaseLib/ConfigTree.h"
-#include "NumLib/ODESolver/NonlinearSolver.h"
-#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
-#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
-
-#include "ProjectData.h"
-
-namespace ApplicationsLib
-{
-//! Time loop capable of time-integrating several uncoupled processes at once.
-class UncoupledProcessesTimeLoop
-{
-public:
-    explicit UncoupledProcessesTimeLoop(
-        std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper)
-        : _timestepper{std::move(timestepper)}
-    {
-    }
-
-    bool loop(ProjectData& project);
-
-    ~UncoupledProcessesTimeLoop();
-
-private:
-    //! An abstract nonlinear solver
-    using AbstractNLSolver = NumLib::NonlinearSolverBase;
-    //! An abstract equations system
-    using EquationSystem = NumLib::EquationSystem;
-    //! An abstract process
-    using Process = ProcessLib::Process;
-    //! An abstract time discretization
-    using TimeDisc = NumLib::TimeDiscretization;
-
-    std::vector<GlobalVector*> _process_solutions;
-    std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
-
-    struct SingleProcessData
-    {
-        template <NumLib::ODESystemTag ODETag, NumLib::NonlinearSolverTag NLTag>
-        SingleProcessData(NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
-                          TimeDisc& time_disc,
-                          NumLib::ODESystem<ODETag, NLTag>& ode_sys)
-            : nonlinear_solver_tag(NLTag),
-              nonlinear_solver(nonlinear_solver),
-              tdisc_ode_sys(new NumLib::TimeDiscretizedODESystem<ODETag, NLTag>(
-                  ode_sys, time_disc)),
-              mat_strg(
-                  dynamic_cast<NumLib::InternalMatrixStorage&>(*tdisc_ode_sys))
-        {
-        }
-
-        SingleProcessData(SingleProcessData&& spd)
-            : nonlinear_solver_tag(spd.nonlinear_solver_tag),
-              nonlinear_solver(spd.nonlinear_solver),
-              tdisc_ode_sys(std::move(spd.tdisc_ode_sys)),
-              mat_strg(spd.mat_strg)
-        {
-        }
-
-        //! Tag containing the missing type information necessary to cast the
-        //! other members of this struct to their concrety types.
-        NumLib::NonlinearSolverTag const nonlinear_solver_tag;
-        AbstractNLSolver& nonlinear_solver;
-        //! type-erased time-discretized ODE system
-        std::unique_ptr<EquationSystem> tdisc_ode_sys;
-        //! cast of \c tdisc_ode_sys to NumLib::InternalMatrixStorage
-        NumLib::InternalMatrixStorage& mat_strg;
-    };
-
-    /*! Creates a new instance of PerProcessData from the given data.
-     *
-     * \tparam ODETag tag indicating the type of ODE to be solved.
-     *
-     * \param nonlinear_solver the nonlinear solver to be used
-     * \param ode_sys   the ODE (i.e. the process) to be integrated in time
-     * \param time_disc the time discretization used for integrating \c ode_sys
-     *
-     * \note Currently the \c ODETag is automatically inferred from the given
-     *       \c ody_sys. This works as long as \c Process derives from
-     *       \c ODESystem<GlobalMatrix, GlobalVector, ODETag>, i.e. as long we
-     *       only deal with one type of ODE. When we introduce more types, this
-     *       method will have to be extended slightly.
-     */
-    template <NumLib::ODESystemTag ODETag>
-    SingleProcessData makeSingleProcessData(
-        AbstractNLSolver& nonlinear_solver,
-        NumLib::ODESystem<ODETag, NumLib::NonlinearSolverTag::Picard>& ode_sys,
-        TimeDisc& time_disc)
-    {
-        using Tag = NumLib::NonlinearSolverTag;
-        // A concrete Picard solver
-        using NonlinearSolverPicard = NumLib::NonlinearSolver<Tag::Picard>;
-        // A concrete Newton solver
-        using NonlinearSolverNewton = NumLib::NonlinearSolver<Tag::Newton>;
-
-        if (auto* nonlinear_solver_picard =
-                dynamic_cast<NonlinearSolverPicard*>(&nonlinear_solver))
-        {
-            // The Picard solver can also work with a Newton-ready ODE,
-            // because the Newton ODESystem derives from the Picard ODESystem.
-            // So no further checks are needed here.
-
-            return SingleProcessData{*nonlinear_solver_picard, time_disc,
-                                     ode_sys};
-        }
-        else if (auto* nonlinear_solver_newton =
-                     dynamic_cast<NonlinearSolverNewton*>(&nonlinear_solver))
-        {
-            // The Newton-Raphson method needs a Newton-ready ODE.
-
-            using ODENewton = NumLib::ODESystem<ODETag, Tag::Newton>;
-            if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
-            {
-                return SingleProcessData{*nonlinear_solver_newton, time_disc,
-                                         *ode_newton};
-            }
-            else
-            {
-                OGS_FATAL(
-                    "You are trying to solve a non-Newton-ready ODE with the"
-                    " Newton-Raphson method. Aborting");
-            }
-        }
-        else
-        {
-            OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
-        }
-    }
-
-    //! Constructs and returns a \c vector of SingleProcessData from the given
-    //! \c project.
-    std::vector<SingleProcessData> initInternalData(ProjectData& project);
-
-    //! Sets initial conditions for the given \c project and \c
-    //! per_process_data.
-    void setInitialConditions(ProjectData& project, double const t0,
-                              std::vector<SingleProcessData>& per_process_data);
-
-    //! Solves one timestep for the given \c process.
-    bool solveOneTimeStepOneProcess(
-            GlobalVector& x, std::size_t const timestep, double const t, double const delta_t,
-            SingleProcessData& process_data,
-            Process& process, ProcessLib::Output const& output_control);
-
-    //! Sets the EquationSystem for the given nonlinear solver,
-    //! which is Picard or Newton depending on the NLTag.
-    template <NumLib::NonlinearSolverTag NLTag>
-    static void setEquationSystem(AbstractNLSolver& nonlinear_solver,
-                                  EquationSystem& eq_sys,
-                                  NumLib::ConvergenceCriterion& conv_crit)
-    {
-        using Solver = NumLib::NonlinearSolver<NLTag>;
-        using EqSys = NumLib::NonlinearSystem<NLTag>;
-
-        assert(dynamic_cast<Solver*>(&nonlinear_solver) != nullptr);
-        assert(dynamic_cast<EqSys*>(&eq_sys) != nullptr);
-
-        auto& nl_solver_ = static_cast<Solver&>(nonlinear_solver);
-        auto& eq_sys_ = static_cast<EqSys&>(eq_sys);
-
-        nl_solver_.setEquationSystem(eq_sys_, conv_crit);
-    }
-
-    //! Sets the EquationSystem for the given nonlinear solver,
-    //! transparently both for Picard and Newton solvers.
-    static void setEquationSystem(AbstractNLSolver& nonlinear_solver,
-                                  EquationSystem& eq_sys,
-                                  NumLib::ConvergenceCriterion& conv_crit,
-                                  NumLib::NonlinearSolverTag nl_tag)
-    {
-        using Tag = NumLib::NonlinearSolverTag;
-        switch (nl_tag)
-        {
-            case Tag::Picard:
-                setEquationSystem<Tag::Picard>(nonlinear_solver, eq_sys,
-                                               conv_crit);
-                break;
-            case Tag::Newton:
-                setEquationSystem<Tag::Newton>(nonlinear_solver, eq_sys,
-                                               conv_crit);
-                break;
-        }
-    }
-
-    //! Applies known solutions to the solution vector \c x, transparently
-    //! for equation systems linearized with either the Picard or Newton method.
-    template <NumLib::NonlinearSolverTag NLTag>
-    static void applyKnownSolutions(EquationSystem const& eq_sys,
-                                    GlobalVector& x)
-    {
-        using EqSys = NumLib::NonlinearSystem<NLTag>;
-        assert(dynamic_cast<EqSys const*>(&eq_sys) != nullptr);
-        auto& eq_sys_ = static_cast<EqSys const&>(eq_sys);
-
-        eq_sys_.applyKnownSolutions(x);
-    }
-
-    //! Applies known solutions to the solution vector \c x, transparently
-    //! for equation systems linearized with either the Picard or Newton method.
-    static void applyKnownSolutions(EquationSystem const& eq_sys,
-                                    NumLib::NonlinearSolverTag const nl_tag,
-                                    GlobalVector& x)
-    {
-        using Tag = NumLib::NonlinearSolverTag;
-        switch (nl_tag)
-        {
-            case Tag::Picard:
-                applyKnownSolutions<Tag::Picard>(eq_sys, x);
-                break;
-            case Tag::Newton:
-                applyKnownSolutions<Tag::Newton>(eq_sys, x);
-                break;
-        }
-    }
-};
-
-//! Builds an UncoupledProcessesTimeLoop from the given configuration.
-std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
-    BaseLib::ConfigTree const& conf);
-
-}  // namespace ApplicationsLib
-
-#endif  // APPLICATIONSLIB_UNCOUPLED_PROCESSES_TIMELOOP
diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index 3fb12e57e45c59b5f9fc30e4d6d2c7753a1f3808..e6a9087ba2c892e894515bfe90e159d193fc25c7 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -23,7 +23,7 @@
 #include "Applications/ApplicationsLib/LinearSolverLibrarySetup.h"
 #include "Applications/ApplicationsLib/LogogSetup.h"
 #include "Applications/ApplicationsLib/ProjectData.h"
-#include "Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h"
+#include "ProcessLib/UncoupledProcessesTimeLoop.h"
 
 #include "NumLib/NumericsConfig.h"
 
@@ -108,16 +108,12 @@ int main(int argc, char *argv[])
             project_config.checkAndInvalidate();
             BaseLib::ConfigTree::assertNoSwallowedErrors();
 
-            // Create processes.
-            project.buildProcesses();
-
             BaseLib::ConfigTree::assertNoSwallowedErrors();
 
             INFO("Initialize processes.");
-            for (auto p_it = project.processesBegin();
-                 p_it != project.processesEnd(); ++p_it)
+            for (auto& p : project.getProcesses())
             {
-                (*p_it)->initialize();
+                p.second->initialize();
             }
 
             BaseLib::ConfigTree::assertNoSwallowedErrors();
@@ -125,7 +121,7 @@ int main(int argc, char *argv[])
             INFO("Solve processes.");
 
             auto& time_loop = project.getTimeLoop();
-            solver_succeeded = time_loop.loop(project);
+            solver_succeeded = time_loop.loop();
         }  // This nested scope ensures that everything that could possibly
            // possess a ConfigTree is destructed before the final check below is
            // done.
diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
index 64b79443c5e386ccdd158186fff3b74492317777..f24d8e78510a2a88156f62e9133df2341a8d023b 100644
--- a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
@@ -20,9 +20,6 @@ namespace GroundwaterFlow
 {
 std::unique_ptr<Process> createGroundwaterFlowProcess(
     MeshLib::Mesh& mesh,
-    Process::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config)
@@ -58,15 +55,9 @@ std::unique_ptr<Process> createGroundwaterFlowProcess(
     ProcessLib::parseSecondaryVariables(config, secondary_variables,
                                         named_function_caller);
 
-    //! \ogs_file_param{process__output}
-    ProcessOutput process_output{config.getConfigSubtree("output")};
-
     return std::unique_ptr<Process>{new GroundwaterFlowProcess{
-        mesh, nonlinear_solver, std::move(time_discretization),
-        std::move(convergence_criterion), parameters,
-        std::move(process_variables), std::move(process_data),
-        std::move(secondary_variables), std::move(process_output),
-        std::move(named_function_caller)}};
+        mesh, parameters, std::move(process_variables), std::move(process_data),
+        std::move(secondary_variables), std::move(named_function_caller)}};
 }
 
 }  // namespace GroundwaterFlow
diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h
index 3e964c246f8653e4963fc56ed828e7f01e80108e..0f3c11ce5116147522410d544b90fc97ad4016fe 100644
--- a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h
@@ -20,9 +20,6 @@ namespace GroundwaterFlow
 {
 std::unique_ptr<Process> createGroundwaterFlowProcess(
     MeshLib::Mesh& mesh,
-    Process::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index e077af278382462a8f329036176aa8299f5a19d9..fd13527d1eacff4394ee6da8a6ea3044bf5bfc55 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -19,33 +19,15 @@ namespace GroundwaterFlow
 {
 GroundwaterFlowProcess::GroundwaterFlowProcess(
     MeshLib::Mesh& mesh,
-    Base::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Base::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     GroundwaterFlowProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
-    ProcessOutput&& process_output,
     NumLib::NamedFunctionCaller&& named_function_caller)
-    : Process(mesh, nonlinear_solver, std::move(time_discretization),
-              std::move(convergence_criterion), parameters,
-              std::move(process_variables), std::move(secondary_variables),
-              std::move(process_output), std::move(named_function_caller)),
+    : Process(mesh, parameters, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller)),
       _process_data(std::move(process_data))
 {
-    if (dynamic_cast<NumLib::ForwardEuler*>(
-            &Base::getTimeDiscretization()) != nullptr)
-    {
-        OGS_FATAL(
-            "GroundwaterFlowProcess can not be solved with the ForwardEuler"
-            " time discretization scheme. Aborting");
-        // Because the M matrix is not assembled. Thus, the linearized system
-        // would be singular. The same applies to CrankNicolson with theta = 0.0,
-        // but this case is not checked here.
-        // Anyway, the GroundwaterFlowProcess shall be transferred to a simpler
-        // ODESystemTag in the future.
-    }
 }
 
 void GroundwaterFlowProcess::initializeConcreteProcess(
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 8f043ea6dcde6275957a9bb7be7f4a23be563afd..ca21d362461104aee9c49d2d0a32782286ae8cb5 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -27,15 +27,11 @@ class GroundwaterFlowProcess final : public Process
 public:
     GroundwaterFlowProcess(
         MeshLib::Mesh& mesh,
-        Base::NonlinearSolver& nonlinear_solver,
-        std::unique_ptr<Base::TimeDiscretization>&& time_discretization,
-        std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         std::vector<std::reference_wrapper<ProcessVariable>>&&
             process_variables,
         GroundwaterFlowProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,
-        ProcessOutput&& process_output,
         NumLib::NamedFunctionCaller&& named_function_caller);
 
     //! \name ODESystem interface
diff --git a/ProcessLib/Output.cpp b/ProcessLib/Output.cpp
index fa94ca94ce7898da61b0a7a0456fb6584c8e3ad2..24c871b8f6a78fb7c1ff2441fc8e9fbd71fbf56c 100644
--- a/ProcessLib/Output.cpp
+++ b/ProcessLib/Output.cpp
@@ -86,25 +86,20 @@ newInstance(const BaseLib::ConfigTree &config, std::string const& output_directo
     return out;
 }
 
-void Output::
-initialize(Output::ProcessIter first, const Output::ProcessIter& last)
+void Output::addProcess(ProcessLib::Process const& process, const unsigned pcs_idx)
 {
-    for (unsigned pcs_idx = 0; first != last; ++first, ++pcs_idx)
-    {
-        auto const filename = _output_file_prefix
-                              + "_pcs_" + std::to_string(pcs_idx)
-                              + ".pvd";
-        _single_process_data.emplace(std::piecewise_construct,
-                std::forward_as_tuple(&**first),
-                std::forward_as_tuple(pcs_idx, filename));
-    }
+    auto const filename =
+        _output_file_prefix + "_pcs_" + std::to_string(pcs_idx) + ".pvd";
+    _single_process_data.emplace(std::piecewise_construct,
+                                 std::forward_as_tuple(&process),
+                                 std::forward_as_tuple(pcs_idx, filename));
 }
 
-void Output::
-doOutputAlways(Process const& process,
-               unsigned timestep,
-               const double t,
-               GlobalVector const& x)
+void Output::doOutputAlways(Process const& process,
+                            ProcessOutput const& process_output,
+                            unsigned timestep,
+                            const double t,
+                            GlobalVector const& x)
 {
     BaseLib::RunTime time_output;
     time_output.start();
@@ -122,33 +117,36 @@ doOutputAlways(Process const& process,
             + "_t_"  + std::to_string(t)
             + ".vtu";
     DBUG("output to %s", output_file_name.c_str());
-    process.output(output_file_name, timestep, x);
+    doProcessOutput(output_file_name, x, process.getMesh(),
+                    process.getDOFTable(), process.getProcessVariables(),
+                    process.getSecondaryVariables(), process_output);
     spd.pvd_file.addVTUFile(output_file_name, t);
 
     INFO("[time] Output took %g s.", time_output.elapsed());
 }
 
-void Output::
-doOutput(Process const& process,
-         unsigned timestep,
-         const double t,
-         GlobalVector const& x)
+void Output::doOutput(Process const& process,
+                      ProcessOutput const& process_output,
+                      unsigned timestep,
+                      const double t,
+                      GlobalVector const& x)
 {
     if (shallDoOutput(timestep, _repeats_each_steps))
-        doOutputAlways(process, timestep, t, x);
+        doOutputAlways(process, process_output, timestep, t, x);
 }
 
-void Output::
-doOutputLastTimestep(Process const& process,
-                     unsigned timestep,
-                     const double t,
-                     GlobalVector const& x)
+void Output::doOutputLastTimestep(Process const& process,
+                                  ProcessOutput const& process_output,
+                                  unsigned timestep,
+                                  const double t,
+                                  GlobalVector const& x)
 {
     if (!shallDoOutput(timestep, _repeats_each_steps))
-        doOutputAlways(process, timestep, t, x);
+        doOutputAlways(process, process_output, timestep, t, x);
 }
 
 void Output::doOutputNonlinearIteration(Process const& process,
+                                        ProcessOutput const& process_output,
                                         const unsigned timestep, const double t,
                                         GlobalVector const& x,
                                         const unsigned iteration) const
@@ -172,7 +170,9 @@ void Output::doOutputNonlinearIteration(Process const& process,
             + "_nliter_" + std::to_string(iteration)
             + ".vtu";
     DBUG("output iteration results to %s", output_file_name.c_str());
-    process.output(output_file_name, timestep, x);
+    doProcessOutput(output_file_name, x, process.getMesh(),
+                    process.getDOFTable(), process.getProcessVariables(),
+                    process.getSecondaryVariables(), process_output);
 
     INFO("[time] Output took %g s.", time_output.elapsed());
 }
diff --git a/ProcessLib/Output.h b/ProcessLib/Output.h
index a5284378f9eadf825b937733c5f0e6ed1e7baa72..8071b2cfadbdad0f1802c89b6d3cbdd2cccb01e5 100644
--- a/ProcessLib/Output.h
+++ b/ProcessLib/Output.h
@@ -13,6 +13,7 @@
 #include "BaseLib/ConfigTree.h"
 #include "MeshLib/IO/VtkIO/PVDFile.h"
 #include "Process.h"
+#include "ProcessOutput.h"
 
 namespace ProcessLib
 {
@@ -29,38 +30,33 @@ public:
     newInstance(const BaseLib::ConfigTree& config,
                 const std::string& output_directory);
 
-    using ProcessIter = std::vector<std::unique_ptr<ProcessLib::Process>>
-                        ::const_iterator;
-
-    //! Opens a PVD file for each process.
-    void initialize(ProcessIter first, const ProcessIter& last);
+    //! TODO doc. Opens a PVD file for each process.
+    void addProcess(ProcessLib::Process const& process, const unsigned pcs_idx);
 
     //! Writes output for the given \c process if it should be written in the
     //! given \c timestep.
-    void doOutput(
-            Process const& process, unsigned timestep,
-            const double t,
-            GlobalVector const& x);
+    void doOutput(Process const& process, ProcessOutput const& process_output,
+                  unsigned timestep, const double t, GlobalVector const& x);
 
     //! Writes output for the given \c process if it has not been written yet.
     //! This method is intended for doing output after the last timestep in order
     //! to make sure that its results are written.
-    void doOutputLastTimestep(
-            Process const& process, unsigned timestep,
-            const double t,
-            GlobalVector const& x);
+    void doOutputLastTimestep(Process const& process,
+                              ProcessOutput const& process_output,
+                              unsigned timestep, const double t,
+                              GlobalVector const& x);
 
     //! Writes output for the given \c process.
     //! This method will always write.
     //! It is intended to write output in error handling routines.
-    void doOutputAlways(
-            Process const& process, unsigned timestep,
-            const double t,
-            GlobalVector const& x);
+    void doOutputAlways(Process const& process,
+                        ProcessOutput const& process_output, unsigned timestep,
+                        const double t, GlobalVector const& x);
 
     //! Writes output for the given \c process.
     //! To be used for debug output after an iteration of the nonlinear solver.
     void doOutputNonlinearIteration(Process const& process,
+                                    ProcessOutput const& process_output,
                                     const unsigned timestep, const double t,
                                     GlobalVector const& x,
                                     const unsigned iteration) const;
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 301edeb82ac254e0f3090dbf2012e056eff11c6a..1b325595532a864394f2159c87852dc44c89eee1 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -20,34 +20,18 @@ namespace ProcessLib
 {
 Process::Process(
     MeshLib::Mesh& mesh,
-    NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     SecondaryVariableCollection&& secondary_variables,
-    ProcessOutput&& process_output,
     NumLib::NamedFunctionCaller&& named_function_caller)
     : _mesh(mesh),
       _secondary_variables(std::move(secondary_variables)),
-      _process_output(std::move(process_output)),
       _named_function_caller(std::move(named_function_caller)),
-      _nonlinear_solver(nonlinear_solver),
-      _time_discretization(std::move(time_discretization)),
-      _convergence_criterion(std::move(convergence_criterion)),
       _process_variables(std::move(process_variables)),
       _boundary_conditions(parameters)
 {
 }
 
-void Process::output(std::string const& file_name,
-                     const unsigned /*timestep*/,
-                     GlobalVector const& x) const
-{
-    doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map,
-                    _process_variables, _secondary_variables, _process_output);
-}
-
 void Process::initialize()
 {
     DBUG("Initialize process.");
@@ -192,11 +176,6 @@ void Process::constructDofTable()
     _local_to_global_index_map.reset(new NumLib::LocalToGlobalIndexMap(
         std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION));
 
-    if (auto* conv_crit =
-            dynamic_cast<NumLib::ConvergenceCriterionPerComponent*>(
-                _convergence_criterion.get())) {
-        conv_crit->setDOFTable(*_local_to_global_index_map, _mesh);
-    }
 }
 
 void Process::initializeExtrapolator()
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 650ffb7d6fbb6067f2aaa08f592a0d52981d7cf1..d37975c072eb0c58ee860b45132d9c7c27d6d1ed 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -18,7 +18,7 @@
 #include "ProcessLib/Parameter/Parameter.h"
 
 #include "ExtrapolatorData.h"
-#include "ProcessOutput.h"
+#include "ProcessVariable.h"
 #include "SecondaryVariable.h"
 #include "CachedSecondaryVariable.h"
 
@@ -40,14 +40,10 @@ public:
 
     Process(
         MeshLib::Mesh& mesh,
-        NonlinearSolver& nonlinear_solver,
-        std::unique_ptr<TimeDiscretization>&& time_discretization,
-        std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         std::vector<std::reference_wrapper<ProcessVariable>>&&
             process_variables,
         SecondaryVariableCollection&& secondary_variables,
-        ProcessOutput&& process_output,
         NumLib::NamedFunctionCaller&& named_function_caller);
 
     /// Preprocessing before starting assembly for new timestep.
@@ -64,12 +60,6 @@ public:
 
     NumLib::IterationResult postIteration(GlobalVector const& x) override final;
 
-    /// Process output.
-    /// The file_name is indicating the name of possible output file.
-    void output(std::string const& file_name,
-                const unsigned /*timestep*/,
-                GlobalVector const& x) const;
-
     void initialize();
 
     void setInitialConditions(const double t, GlobalVector& x);
@@ -93,15 +83,22 @@ public:
         return _boundary_conditions.getKnownSolutions(t);
     }
 
-    NonlinearSolver& getNonlinearSolver() const { return _nonlinear_solver; }
-    TimeDiscretization& getTimeDiscretization() const
+    NumLib::LocalToGlobalIndexMap const& getDOFTable() const
+    {
+        return *_local_to_global_index_map;
+    }
+
+    MeshLib::Mesh& getMesh() const { return _mesh; }
+
+    std::vector<std::reference_wrapper<ProcessVariable>> const&
+    getProcessVariables() const
     {
-        return *_time_discretization;
+        return _process_variables;
     }
 
-    NumLib::ConvergenceCriterion& getConvergenceCriterion() const
+    SecondaryVariableCollection const& getSecondaryVariables() const
     {
-        return *_convergence_criterion;
+        return _secondary_variables;
     }
 
 protected:
@@ -160,7 +157,6 @@ protected:
     std::unique_ptr<NumLib::LocalToGlobalIndexMap> _local_to_global_index_map;
 
     SecondaryVariableCollection _secondary_variables;
-    ProcessOutput _process_output;
 
     NumLib::NamedFunctionCaller _named_function_caller;
     std::vector<std::unique_ptr<CachedSecondaryVariable>>
@@ -171,10 +167,6 @@ private:
     unsigned const _integration_order = 2;
     GlobalSparsityPattern _sparsity_pattern;
 
-    NonlinearSolver& _nonlinear_solver;
-    std::unique_ptr<TimeDiscretization> _time_discretization;
-    std::unique_ptr<NumLib::ConvergenceCriterion> _convergence_criterion;
-
     /// Variables used by this process.
     std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
 
diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp
index a1b6c90322e8b364e756a9d5a551b1a999aaa92e..f57906d55806ae39537fccc2dcbba194eb6296b9 100644
--- a/ProcessLib/ProcessOutput.cpp
+++ b/ProcessLib/ProcessOutput.cpp
@@ -38,16 +38,6 @@ ProcessOutput::ProcessOutput(BaseLib::ConfigTree const& output_config)
     {
         output_residuals = *out_resid;
     }
-
-    // debug output
-    if (auto const param =
-        //! \ogs_file_param{process__output__output_iteration_results}
-        output_config.getConfigParameterOptional<bool>("output_iteration_results"))
-    {
-        DBUG("output_iteration_results: %s", (*param) ? "true" : "false");
-
-        output_iteration_results = *param;
-    }
 }
 
 
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
index db13504f76676e20efd38b0c4abb6655586bd8fe..44b40f4459ea8dbead96e6fdb778baf674e68dfd 100644
--- a/ProcessLib/ProcessOutput.h
+++ b/ProcessLib/ProcessOutput.h
@@ -27,8 +27,6 @@ struct ProcessOutput final
 
     //! Tells if also to output extrapolation residuals.
     bool output_residuals = false;
-
-    bool output_iteration_results = false;
 };
 
 
diff --git a/ProcessLib/TES/CreateTESProcess.cpp b/ProcessLib/TES/CreateTESProcess.cpp
index 2c5393eb0d673429a42fac73603b236a2b7ca5c1..0f2b0da925436ea2f9c065cea4810519ace9a775 100644
--- a/ProcessLib/TES/CreateTESProcess.cpp
+++ b/ProcessLib/TES/CreateTESProcess.cpp
@@ -18,9 +18,6 @@ namespace TES
 {
 std::unique_ptr<Process> createTESProcess(
     MeshLib::Mesh& mesh,
-    Process::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     BaseLib::ConfigTree const& config)
@@ -42,14 +39,9 @@ std::unique_ptr<Process> createTESProcess(
     ProcessLib::parseSecondaryVariables(config, secondary_variables,
                                         named_function_caller);
 
-    //! \ogs_file_param{process__output}
-    ProcessOutput process_output{config.getConfigSubtree("output")};
-
     return std::unique_ptr<Process>{new TESProcess{
-        mesh, nonlinear_solver, std::move(time_discretization),
-        std::move(convergence_criterion), parameters,
-        std::move(process_variables), std::move(secondary_variables),
-        std::move(process_output), std::move(named_function_caller), config}};
+        mesh, parameters, std::move(process_variables), std::move(secondary_variables),
+        std::move(named_function_caller), config}};
 }
 
 }  // namespace TES
diff --git a/ProcessLib/TES/CreateTESProcess.h b/ProcessLib/TES/CreateTESProcess.h
index f51d00030e62f638f5625092cccca88c82b9304e..856625928eb0b4e580b30313b5672de13af20019 100644
--- a/ProcessLib/TES/CreateTESProcess.h
+++ b/ProcessLib/TES/CreateTESProcess.h
@@ -19,9 +19,6 @@ namespace TES
 {
 std::unique_ptr<Process> createTESProcess(
     MeshLib::Mesh& mesh,
-    Process::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& /*parameters*/,
     BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index c53b16cebdf74d4f1a0bb12cb7e95d510f03ec6a..40a095841ff8933fb272e50360d3976087a21520 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -55,19 +55,13 @@ namespace TES
 {
 TESProcess::TESProcess(
     MeshLib::Mesh& mesh,
-    Process::NonlinearSolver& nonlinear_solver,
-    std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-    std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     SecondaryVariableCollection&& secondary_variables,
-    ProcessOutput&& process_output,
     NumLib::NamedFunctionCaller&& named_function_caller,
     const BaseLib::ConfigTree& config)
-    : Process(mesh, nonlinear_solver, std::move(time_discretization),
-              std::move(convergence_criterion), parameters,
-              std::move(process_variables), std::move(secondary_variables),
-              std::move(process_output), std::move(named_function_caller))
+    : Process(mesh, parameters, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller))
 {
     DBUG("Create TESProcess.");
 
@@ -298,20 +292,6 @@ void TESProcess::preIterationConcreteProcess(const unsigned iter,
 NumLib::IterationResult TESProcess::postIterationConcreteProcess(
     GlobalVector const& x)
 {
-    if (this->_process_output.output_iteration_results)
-    {
-        DBUG("output results of iteration %li",
-             _assembly_params.total_iteration);
-        std::string fn =
-            "tes_iter_" + std::to_string(_assembly_params.total_iteration) +
-            +"_ts_" + std::to_string(_assembly_params.timestep) + "_" +
-            std::to_string(_assembly_params.iteration_in_current_timestep) +
-            "_" + std::to_string(_assembly_params.number_of_try_of_iteration) +
-            ".vtu";
-
-        this->output(fn, 0, x);
-    }
-
     bool check_passed = true;
 
     if (!Trafo::constrained)
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 84dbb50132ee213c17810cb1c3eae68fcfd2f4af..e0d94642f6f36b0861c52ee15714afb67962d8f8 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -33,14 +33,10 @@ class TESProcess final : public Process
 public:
     TESProcess(
         MeshLib::Mesh& mesh,
-        Process::NonlinearSolver& nonlinear_solver,
-        std::unique_ptr<Process::TimeDiscretization>&& time_discretization,
-        std::unique_ptr<NumLib::ConvergenceCriterion>&& convergence_criterion,
         std::vector<std::unique_ptr<ParameterBase>> const& parameters,
         std::vector<std::reference_wrapper<ProcessVariable>>&&
             process_variables,
         SecondaryVariableCollection&& secondary_variables,
-        ProcessOutput&& process_output,
         NumLib::NamedFunctionCaller&& named_function_caller,
         BaseLib::ConfigTree const& config);
 
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4324562487e93638e79968d3e1e167fa8f589066
--- /dev/null
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -0,0 +1,561 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "UncoupledProcessesTimeLoop.h"
+#include "BaseLib/uniqueInsert.h"
+#include "BaseLib/RunTime.h"
+#include "NumLib/ODESolver/TimeDiscretizationBuilder.h"
+#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
+#include "NumLib/ODESolver/ConvergenceCriterionPerComponent.h"
+#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+
+std::unique_ptr<NumLib::ITimeStepAlgorithm> createTimeStepper(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__time_loop__time_stepping__type}
+    auto const type = config.peekConfigParameter<std::string>("type");
+
+    std::unique_ptr<NumLib::ITimeStepAlgorithm> timestepper;
+
+    if (type == "SingleStep")
+    {
+        //! \ogs_file_param_special{prj__time_stepping__SingleStep}
+        config.ignoreConfigParameter("type");
+        timestepper.reset(new NumLib::FixedTimeStepping(0.0, 1.0, 1.0));
+    }
+    else if (type == "FixedTimeStepping")
+    {
+        timestepper = NumLib::FixedTimeStepping::newInstance(config);
+    }
+    else
+    {
+        OGS_FATAL("Unknown timestepper type: `%s'.", type.c_str());
+    }
+
+    return timestepper;
+}
+
+std::unique_ptr<ProcessLib::Output> createOutput(
+    BaseLib::ConfigTree const& config, std::string const& output_directory)
+{
+    //! \ogs_file_param{prj__time_loop__output__type}
+    config.checkConfigParameter("type", "VTK");
+    DBUG("Parse output configuration:");
+
+    return ProcessLib::Output::newInstance(config, output_directory);
+}
+
+
+//! Sets the EquationSystem for the given nonlinear solver,
+//! which is Picard or Newton depending on the NLTag.
+template <NumLib::NonlinearSolverTag NLTag>
+static void setEquationSystem(NumLib::NonlinearSolverBase& nonlinear_solver,
+                              NumLib::EquationSystem& eq_sys,
+                              NumLib::ConvergenceCriterion& conv_crit)
+{
+    using Solver = NumLib::NonlinearSolver<NLTag>;
+    using EqSys = NumLib::NonlinearSystem<NLTag>;
+
+    assert(dynamic_cast<Solver*>(&nonlinear_solver) != nullptr);
+    assert(dynamic_cast<EqSys*>(&eq_sys) != nullptr);
+
+    auto& nl_solver_ = static_cast<Solver&>(nonlinear_solver);
+    auto& eq_sys_ = static_cast<EqSys&>(eq_sys);
+
+    nl_solver_.setEquationSystem(eq_sys_, conv_crit);
+}
+
+//! Sets the EquationSystem for the given nonlinear solver,
+//! transparently both for Picard and Newton solvers.
+static void setEquationSystem(NumLib::NonlinearSolverBase& nonlinear_solver,
+                              NumLib::EquationSystem& eq_sys,
+                              NumLib::ConvergenceCriterion& conv_crit,
+                              NumLib::NonlinearSolverTag nl_tag)
+{
+    using Tag = NumLib::NonlinearSolverTag;
+    switch (nl_tag)
+    {
+        case Tag::Picard:
+            setEquationSystem<Tag::Picard>(nonlinear_solver, eq_sys,
+                                           conv_crit);
+            break;
+        case Tag::Newton:
+            setEquationSystem<Tag::Newton>(nonlinear_solver, eq_sys,
+                                           conv_crit);
+            break;
+    }
+}
+
+//! Applies known solutions to the solution vector \c x, transparently
+//! for equation systems linearized with either the Picard or Newton method.
+template <NumLib::NonlinearSolverTag NLTag>
+static void applyKnownSolutions(NumLib::EquationSystem const& eq_sys,
+                                GlobalVector& x)
+{
+    using EqSys = NumLib::NonlinearSystem<NLTag>;
+    assert(dynamic_cast<EqSys const*>(&eq_sys) != nullptr);
+    auto& eq_sys_ = static_cast<EqSys const&>(eq_sys);
+
+    eq_sys_.applyKnownSolutions(x);
+}
+
+//! Applies known solutions to the solution vector \c x, transparently
+//! for equation systems linearized with either the Picard or Newton method.
+static void applyKnownSolutions(NumLib::EquationSystem const& eq_sys,
+                                NumLib::NonlinearSolverTag const nl_tag,
+                                GlobalVector& x)
+{
+    using Tag = NumLib::NonlinearSolverTag;
+    switch (nl_tag)
+    {
+        case Tag::Picard:
+            applyKnownSolutions<Tag::Picard>(eq_sys, x);
+            break;
+        case Tag::Newton:
+            applyKnownSolutions<Tag::Newton>(eq_sys, x);
+            break;
+    }
+}
+
+namespace ProcessLib
+{
+struct SingleProcessData
+{
+    template <NumLib::NonlinearSolverTag NLTag>
+    SingleProcessData(
+        NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
+        std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
+        std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
+        Process& process_,
+        ProcessOutput&& process_output_);
+
+    SingleProcessData(SingleProcessData&& spd);
+
+    //! Tag containing the missing type information necessary to cast the
+    //! other members of this struct to their concrety types.
+    NumLib::NonlinearSolverTag const nonlinear_solver_tag;
+    NumLib::NonlinearSolverBase& nonlinear_solver;
+    std::unique_ptr<NumLib::ConvergenceCriterion> conv_crit;
+
+    std::unique_ptr<NumLib::TimeDiscretization> time_disc;
+    //! type-erased time-discretized ODE system
+    std::unique_ptr<NumLib::EquationSystem> tdisc_ode_sys;
+    //! cast of \c tdisc_ode_sys to NumLib::InternalMatrixStorage
+    NumLib::InternalMatrixStorage* mat_strg = nullptr;
+
+    Process& process;
+    ProcessOutput process_output;
+};
+
+template <NumLib::NonlinearSolverTag NLTag>
+SingleProcessData::SingleProcessData(
+    NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
+    std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
+    std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
+    Process& process_,
+    ProcessOutput&& process_output_)
+    : nonlinear_solver_tag(NLTag),
+      nonlinear_solver(nonlinear_solver),
+      conv_crit(std::move(conv_crit_)),
+      time_disc(std::move(time_disc_)),
+      process(process_),
+      process_output(std::move(process_output_))
+{
+}
+
+SingleProcessData::SingleProcessData(SingleProcessData&& spd)
+    : nonlinear_solver_tag(spd.nonlinear_solver_tag),
+      nonlinear_solver(spd.nonlinear_solver),
+      conv_crit(std::move(spd.conv_crit)),
+      time_disc(std::move(spd.time_disc)),
+      tdisc_ode_sys(std::move(spd.tdisc_ode_sys)),
+      mat_strg(spd.mat_strg),
+      process(spd.process),
+      process_output(std::move(spd.process_output))
+{
+    spd.mat_strg = nullptr;
+}
+
+template <NumLib::ODESystemTag ODETag>
+void setTimeDiscretizedODESystem(
+    SingleProcessData& spd,
+    NumLib::ODESystem<ODETag, NumLib::NonlinearSolverTag::Picard>& ode_sys)
+{
+    using Tag = NumLib::NonlinearSolverTag;
+    // A concrete Picard solver
+    using NonlinearSolverPicard = NumLib::NonlinearSolver<Tag::Picard>;
+    // A concrete Newton solver
+    using NonlinearSolverNewton = NumLib::NonlinearSolver<Tag::Newton>;
+
+    if (dynamic_cast<NonlinearSolverPicard*>(&spd.nonlinear_solver))
+    {
+        // The Picard solver can also work with a Newton-ready ODE,
+        // because the Newton ODESystem derives from the Picard ODESystem.
+        // So no further checks are needed here.
+
+        spd.tdisc_ode_sys.reset(
+            new NumLib::TimeDiscretizedODESystem<ODETag, Tag::Picard>(
+                ode_sys, *spd.time_disc));
+    }
+    else if (dynamic_cast<NonlinearSolverNewton*>(&spd.nonlinear_solver))
+    {
+        // The Newton-Raphson method needs a Newton-ready ODE.
+
+        using ODENewton = NumLib::ODESystem<ODETag, Tag::Newton>;
+        if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
+        {
+            spd.tdisc_ode_sys.reset(
+                new NumLib::TimeDiscretizedODESystem<ODETag, Tag::Newton>(
+                    *ode_newton, *spd.time_disc));
+        }
+        else
+        {
+            OGS_FATAL(
+                "You are trying to solve a non-Newton-ready ODE with the"
+                " Newton-Raphson method. Aborting");
+        }
+    }
+    else
+    {
+        OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
+    }
+
+    spd.mat_strg =
+        dynamic_cast<NumLib::InternalMatrixStorage*>(spd.tdisc_ode_sys.get());
+}
+
+void setTimeDiscretizedODESystem(SingleProcessData& spd)
+{
+    setTimeDiscretizedODESystem(spd, spd.process);
+}
+
+std::unique_ptr<SingleProcessData> makeSingleProcessData(
+    NumLib::NonlinearSolverBase& nonlinear_solver,
+    Process& process,
+    std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
+    std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
+    ProcessOutput&& process_output)
+{
+    using Tag = NumLib::NonlinearSolverTag;
+
+    if (auto* nonlinear_solver_picard =
+            dynamic_cast<NumLib::NonlinearSolver<Tag::Picard>*>(
+                &nonlinear_solver))
+    {
+        return std::unique_ptr<SingleProcessData>{new SingleProcessData{
+            *nonlinear_solver_picard, std::move(conv_crit),
+            std::move(time_disc), process, std::move(process_output)}};
+    }
+    else if (auto* nonlinear_solver_newton =
+                 dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
+                     &nonlinear_solver))
+    {
+        return std::unique_ptr<SingleProcessData>{new SingleProcessData{
+            *nonlinear_solver_newton, std::move(conv_crit),
+            std::move(time_disc), process, std::move(process_output)}};
+    } else {
+        OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
+    }
+}
+
+std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
+    BaseLib::ConfigTree const& config,
+    const std::map<std::string, std::unique_ptr<Process>>&
+        processes,
+    std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
+        nonlinear_solvers)
+{
+    std::vector<std::unique_ptr<SingleProcessData>> per_process_data;
+
+    //! \ogs_file_param{prj__time_loop__processes__process}
+    for (auto pcs_config : config.getConfigSubtreeList("process"))
+    {
+        //! \ogs_file_attr{prj__time_loop__processes__process__ref}
+        auto const pcs_name = pcs_config.getConfigAttribute<std::string>("ref");
+        auto& pcs = *BaseLib::getOrError(
+            processes, pcs_name,
+            "A process with the given name has not been defined.");
+
+        //! \ogs_file_param{prj__time_loop__processes__process__nonlinear_solver}
+        auto const nl_slv_name =
+            pcs_config.getConfigParameter<std::string>("nonlinear_solver");
+        auto& nl_slv = *BaseLib::getOrError(
+            nonlinear_solvers, nl_slv_name,
+            "A nonlinear solver with the given name has not been defined.");
+
+        auto time_disc = NumLib::createTimeDiscretization(
+            //! \ogs_file_param{prj__time_loop__processes__process__time_discretization}
+            pcs_config.getConfigSubtree("time_discretization"));
+
+        auto conv_crit = NumLib::createConvergenceCriterion(
+            //! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
+            pcs_config.getConfigSubtree("convergence_criterion"));
+
+        //! \ogs_file_param{prj__time_loop__processes__process__output}
+        ProcessOutput process_output{
+            pcs_config.getConfigSubtree("output")};
+
+        per_process_data.emplace_back(makeSingleProcessData(
+            nl_slv, pcs, std::move(time_disc), std::move(conv_crit),
+            std::move(process_output)));
+    }
+
+    if (per_process_data.size() != processes.size())
+        OGS_FATAL(
+            "Some processes have not been configured to be solved by this time "
+            "time loop.");
+
+    return per_process_data;
+}
+
+std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
+    BaseLib::ConfigTree const& config, std::string const& output_directory,
+    const std::map<std::string, std::unique_ptr<Process>>&
+        processes,
+    const std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>>&
+        nonlinear_solvers)
+{
+    //! \ogs_file_param{prj__time_loop__time_stepping}
+    auto timestepper =
+        createTimeStepper(config.getConfigSubtree("time_stepping"));
+
+    //! \ogs_file_param{prj__time_loop__output}
+    auto output =
+        createOutput(config.getConfigSubtree("output"), output_directory);
+
+    //! \ogs_file_param{prj__time_loop__processes}
+    auto per_process_data = createPerProcessData(
+        config.getConfigSubtree("processes"), processes, nonlinear_solvers);
+
+    return std::unique_ptr<UncoupledProcessesTimeLoop>{
+        new UncoupledProcessesTimeLoop{std::move(timestepper),
+                                       std::move(output),
+                                       std::move(per_process_data)}};
+}
+
+std::vector<GlobalVector*> setInitialConditions(
+    double const t0,
+    std::vector<std::unique_ptr<SingleProcessData>> const& per_process_data)
+{
+    std::vector<GlobalVector*> process_solutions;
+
+    unsigned pcs_idx = 0;
+    for (auto& spd : per_process_data)
+    {
+        auto& pcs = spd->process;
+        auto& time_disc = *spd->time_disc;
+
+        auto& ode_sys = *spd->tdisc_ode_sys;
+        auto const nl_tag = spd->nonlinear_solver_tag;
+
+        // append a solution vector of suitable size
+        process_solutions.emplace_back(
+            &NumLib::GlobalVectorProvider::provider.getVector(
+                ode_sys.getMatrixSpecifications()));
+
+        auto& x0 = *process_solutions[pcs_idx];
+        pcs.setInitialConditions(t0, x0);
+        MathLib::LinAlg::finalizeAssembly(x0);
+
+        time_disc.setInitialState(t0, x0);  // push IC
+
+        if (time_disc.needsPreload())
+        {
+            auto& nonlinear_solver = spd->nonlinear_solver;
+            auto& mat_strg = *spd->mat_strg;
+            auto& conv_crit = *spd->conv_crit;
+
+            setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
+            nonlinear_solver.assemble(x0);
+            time_disc.pushState(
+                t0, x0, mat_strg);  // TODO: that might do duplicate work
+        }
+
+        ++pcs_idx;
+    }
+
+    return process_solutions;
+}
+
+bool solveOneTimeStepOneProcess(GlobalVector& x, std::size_t const timestep,
+                                double const t, double const delta_t,
+                                SingleProcessData& process_data,
+                                Output const& output_control)
+{
+    auto& process = process_data.process;
+    auto& time_disc = *process_data.time_disc;
+    auto& conv_crit = *process_data.conv_crit;
+    auto& ode_sys = *process_data.tdisc_ode_sys;
+    auto& nonlinear_solver = process_data.nonlinear_solver;
+    auto const nl_tag = process_data.nonlinear_solver_tag;
+
+    setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
+
+    // Note: Order matters!
+    // First advance to the next timestep, then set known solutions at that
+    // time, afterwards pass the right solution vector and time to the
+    // preTimestep() hook.
+
+    time_disc.nextTimestep(t, delta_t);
+
+    applyKnownSolutions(ode_sys, nl_tag, x);
+
+    process.preTimestep(x, t, delta_t);
+
+    auto const post_iteration_callback = [&](
+        unsigned iteration, GlobalVector const& x) {
+        output_control.doOutputNonlinearIteration(
+            process, process_data.process_output, timestep, t, x, iteration);
+    };
+
+    bool nonlinear_solver_succeeded =
+        nonlinear_solver.solve(x, post_iteration_callback);
+
+    auto& mat_strg = *process_data.mat_strg;
+    time_disc.pushState(t, x, mat_strg);
+
+    process.postTimestep(x);
+
+    return nonlinear_solver_succeeded;
+}
+
+UncoupledProcessesTimeLoop::UncoupledProcessesTimeLoop(
+    std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
+    std::unique_ptr<Output>&& output,
+    std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data)
+    : _timestepper{std::move(timestepper)},
+      _output(std::move(output)),
+      _per_process_data(std::move(per_process_data))
+{
+}
+
+bool UncoupledProcessesTimeLoop::loop()
+{
+    // initialize output, convergence criterion, etc.
+    {
+        unsigned pcs_idx = 0;
+        for (auto& spd : _per_process_data)
+        {
+            auto& pcs = spd->process;
+            _output->addProcess(pcs, pcs_idx);
+
+            setTimeDiscretizedODESystem(*spd);
+
+            if (auto* conv_crit =
+                    dynamic_cast<NumLib::ConvergenceCriterionPerComponent*>(
+                        spd->conv_crit.get())) {
+                conv_crit->setDOFTable(pcs.getDOFTable(), pcs.getMesh());
+            }
+
+            ++pcs_idx;
+        }
+    }
+
+    auto const t0 = _timestepper->getTimeStep().current();  // time of the IC
+
+    // init solution storage
+    _process_solutions = setInitialConditions(t0, _per_process_data);
+
+    // output initial conditions
+    {
+        unsigned pcs_idx = 0;
+        for (auto& spd : _per_process_data)
+        {
+            auto& pcs = spd->process;
+            auto const& x0 = *_process_solutions[pcs_idx];
+
+            pcs.preTimestep(x0, t0, _timestepper->getTimeStep().dt());
+            _output->doOutput(pcs, spd->process_output, 0, t0, x0);
+            ++pcs_idx;
+        }
+    }
+
+    double t = t0;
+    std::size_t timestep = 1;  // the first timestep really is number one
+    bool nonlinear_solver_succeeded = true;
+
+    while (_timestepper->next())
+    {
+        BaseLib::RunTime time_timestep;
+        time_timestep.start();
+
+        auto const ts = _timestepper->getTimeStep();
+        auto const delta_t = ts.dt();
+        t = ts.current();
+        timestep = ts.steps();
+
+        INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
+             timestep, t, delta_t);
+
+        // TODO use process name
+        unsigned pcs_idx = 0;
+        for (auto& spd : _per_process_data)
+        {
+            auto& pcs = spd->process;
+            BaseLib::RunTime time_timestep_process;
+            time_timestep_process.start();
+
+            auto& x = *_process_solutions[pcs_idx];
+
+            nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
+                x, timestep, t, delta_t, *spd, *_output);
+
+            INFO("[time] Solving process #%u took %g s in timestep #%u.",
+                 timestep, time_timestep.elapsed(), pcs_idx);
+
+            if (!nonlinear_solver_succeeded)
+            {
+                ERR("The nonlinear solver failed in timestep #%u at t = %g s"
+                    " for process #%u.",
+                    timestep, t, pcs_idx);
+
+                // save unsuccessful solution
+                _output->doOutputAlways(pcs, spd->process_output, timestep, t, x);
+
+                break;
+            }
+            else
+            {
+                _output->doOutput(pcs, spd->process_output, timestep, t, x);
+            }
+
+            ++pcs_idx;
+        }
+
+        INFO("[time] Timestep #%u took %g s.", timestep,
+             time_timestep.elapsed());
+
+        if (!nonlinear_solver_succeeded)
+            break;
+    }
+
+    // output last timestep
+    if (nonlinear_solver_succeeded)
+    {
+        unsigned pcs_idx = 0;
+        for (auto& spd : _per_process_data)
+        {
+            auto& pcs = spd->process;
+            auto const& x = *_process_solutions[pcs_idx];
+            _output->doOutputLastTimestep(pcs, spd->process_output, timestep, t, x);
+
+            ++pcs_idx;
+        }
+    }
+
+    return nonlinear_solver_succeeded;
+}
+
+UncoupledProcessesTimeLoop::~UncoupledProcessesTimeLoop()
+{
+    for (auto* x : _process_solutions)
+        NumLib::GlobalVectorProvider::provider.releaseVector(*x);
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.h b/ProcessLib/UncoupledProcessesTimeLoop.h
new file mode 100644
index 0000000000000000000000000000000000000000..992ab4624e02b1b16c03b20a3a6952e76ce3648d
--- /dev/null
+++ b/ProcessLib/UncoupledProcessesTimeLoop.h
@@ -0,0 +1,56 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_UNCOUPLED_PROCESSES_TIMELOOP
+#define PROCESSLIB_UNCOUPLED_PROCESSES_TIMELOOP
+
+#include <memory>
+#include <logog/include/logog.hpp>
+
+#include "NumLib/ODESolver/NonlinearSolver.h"
+#include "NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h"
+
+#include "Output.h"
+#include "Process.h"
+
+namespace ProcessLib
+{
+struct SingleProcessData;
+
+//! Time loop capable of time-integrating several uncoupled processes at once.
+class UncoupledProcessesTimeLoop
+{
+public:
+    explicit UncoupledProcessesTimeLoop(
+        std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
+        std::unique_ptr<Output>&& output,
+        std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data);
+
+    bool loop();
+
+    ~UncoupledProcessesTimeLoop();
+
+private:
+    std::vector<GlobalVector*> _process_solutions;
+    std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
+    std::unique_ptr<Output> _output;
+    std::vector<std::unique_ptr<SingleProcessData>> _per_process_data;
+};
+
+//! Builds an UncoupledProcessesTimeLoop from the given configuration.
+std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
+    BaseLib::ConfigTree const& config, std::string const& output_directory,
+    std::map<std::string, std::unique_ptr<Process>> const&
+        processes,
+    std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
+        nonlinear_solvers);
+
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_UNCOUPLED_PROCESSES_TIMELOOP
diff --git a/Tests/Data b/Tests/Data
index 08ae549df5b89518a1c8df74bb48ad4e2d583fe1..41e1227f7b04a0ea3e4a520c09a7a2e2c8d3ec28 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 08ae549df5b89518a1c8df74bb48ad4e2d583fe1
+Subproject commit 41e1227f7b04a0ea3e4a520c09a7a2e2c8d3ec28