diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index bfe67f022166e1c23874c7b3a8658c6d81bb3f14..1faa342529f5e560b36115583f664ce207341197 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -5,7 +5,8 @@ The aim of this exercise is to get familiar with setting up coupled free flow/po
 ## Problem set-up
 
 The model domain consists of two non-overlapping subdomains.
-Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium. Both single-phase flow and two-phase flow will be considered in the porous domain.
+Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium.
+Both single-phase flow and two-phase flow will be considered in the porous domain.
 
 
 ### 0. Getting familiar with the code
@@ -14,7 +15,7 @@ Free flow is modeled in the upper subdomain, while the lower subdomain models a
 
 There are three sub folders: `interface` (Exercise 1), `models` (Exercise 2) and `turbulence` (Exercise 3).
 
-The problem-related files for this exercise are:
+The task-related files for the simualtion exercises are as follows:
 * Three __main files__ for the three sub-tasks :`interface/main.cc`, `models/main.cc`, `turbulence/main.cc`,
 * Three __free flow problem files__: `interface/freeflowsubproblem.hh`, `models/freeflowsubproblem.hh`, `turbulence/freeflowsubproblem.hh`
 * Three __porous medium flow problem files__: `interface/porousmediumsubproblem.hh`, `models/porousmediumsubproblem.hh`, `turbulence/porousmediumsubproblem.hh`
@@ -22,26 +23,26 @@ The problem-related files for this exercise are:
 * The __input files__: `interface/params.input`, `models/parmas.input`, `turbulence/params.input`,
 * The __spatial parameters files__: `1pspatialparams.hh`, `2pspatialparams.hh`
 
-In the main file, `TypeTags` for both submodels are defined.
+In the main file, `TypeTags` for both submodels are defined, `FreeflowTypeTag` and `PorousMediumTypeTag`. These `TypeTags` collect all of the properties associated with each subproblem.
 The same applies for types such as `GridManager`, `FVGridGeometry`, `Problem`, etc...
-Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`.
+Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`, which help to assemble and solve the full coupled problem.
 
-The problem files very much look like "regular", uncoupled ones with the exception that they hold a pointer to the `CouplingManager` which allows to evaluate the coupling conditions and to exchange information between the coupled models.
+The problem files very much look like the "regular", uncoupled problem files seen in previous exercises, with the exception that they hold a pointer to the `CouplingManager`.
+This allows them to evaluate the coupling conditions and to exchange information between the coupled models.
 The coupling conditions are realized technically in terms of boundary condition. For instance,
 in `interface/freeflowsubproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow models evaluates the
 mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
 
-Note the certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains
+Note that certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains
 and the sub-control-volume faces at the interface have to match.
 
-
-We will use a staggered grid for the free flow and a cell-centered finite volume method for the porous medium.
+We will use a staggered grid to discretize the free-flow domain and a cell-centered finite volume method for the porous medium domain.
 Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and
-cell center variables (all other variables), therefore in some cases either the `stokesCellCenterIdx`
-or the `stokesFaceIdx` is used respectively, while for the porous medium all variables can be accessed with `darcyIdx`.
+cell center variables (all other variables).
+For this reason either the `CouplingManager::stokesCellCenterIdx` or the `CouplingManager::stokesFaceIdx` index sets need to be used, while for the porous medium all variables can be accessed with `CouplingManager::darcyIdx`.
 
 __Task__:
-Take a closer look at the Stokes/Darcy coupling files before moving to the next part of the exercise:
+Take a closer look at the listed files before moving to the next part of the exercise:
 
 
 ### 1. Changing the interface
@@ -139,7 +140,7 @@ at the position where the coupling boundary conditions are set in `interface/fre
 To check if the simulation behaves as expected, we can compare the velocity profile $`v_x(y)`$ with the analytical solution provided by [Beavers and Joseph (1967)](https://doi.org/10.1017/S0022112067001375).
 For doing so, we uncomment the following line in `main.cc` in the subfolder `interface`.
 ```cpp
-stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
 ```
 
 After re-compiling and re-running the executable, we should be able to see also
@@ -147,23 +148,32 @@ the analytical solution of $`v_x`$ on the free flow domain. Play around with the
 
 __Task C: Change the shape of interface__
 
-Now we want to include a non-flat interface between the two domains. We use `dune-subgrid` to construct two grids for the two domains from one common host grid. To do so, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`. Comment out the first code block and uncomment the second. This will instantiate a host grid and define two helper lambda functions that are used to choose elements from to host grid for the respective sub grid. In the given case, the domain is split in two haves, separated by a sinusoidal interface.
+Now we want to include a non-flat interface between the two domains.
+We use [`dune-subgrid`](https://doi.org/10.1007/s00607-009-0067-2) to construct two grids for the two domains from one common host grid.
+Our hostgrid will be a Dune-Yasp grid (`Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >`)
+and the bounds and resolution of the domain will be set in the `params.input`file under the group `[Grid]`.
+This hostgrid, along with `elementSelector` functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.
+
+To introduce this, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`.
+Comment out the first code block and uncomment the second.
+This will instantiate a host grid and define two helper lambda functions that are used to choose elements from to host grid for the respective sub grid.
+In the given case, the domain is split in two halves, separated by a sinusoidal interface.
 
 ```cpp
-auto elementSelectorStokes = [&](const auto& element)
+auto elementSelectorFreeflow = [&](const auto& element)
 {
-    double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+    double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
     return element.geometry().center()[1] > interface;
 };
 
-auto elementSelectorDarcy = [&](const auto& element)
+auto elementSelectorPorousMedium = [&](const auto& element)
 {
-    double interface  =  params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+    double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
     return element.geometry().center()[1] < interface;
 };
 ```
 
-Make sure, that you have uncommented the lines including the grid managers in both properties files
+Make sure that you have uncommented the lines including the grid managers in both properties files
 ```cpp
 #include <dumux/io/grid/gridmanager_sub.hh>
 ```
@@ -194,11 +204,15 @@ The final result should look something like this:
 
 ![](../extradoc/ex_ff-pm-wave-interface.png)
 
+*Extra Points:*
+Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set.
+What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary?
+Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the `initialAtPos` method.
 
 ### 2. Changing the porous medium model
 
 In this part of the exercise, the coupled system will be extended such that the transport of components
-in both parts is included and the presence of a second phase in the porous medium is considered.
+in both domains is included and the presence of a second phase in the porous medium domain is considered.
 This enables the simulation of the drying of a porous medium (however no energy balance is included yet).
 
 This part of the exercise consists of the following steps:
@@ -230,7 +244,7 @@ liquid saturation as primary variables (p_g-S_l -> `p1s0`) or vice versa.
 
 ```
 template<class TypeTag>
-struct Formulation<TypeTag, TTag::DarcyOnePNC>
+struct Formulation<TypeTag, TTag::PorousMediumOnePNC>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 ```
   in the properties file.
@@ -246,30 +260,41 @@ In this case, the chosen formulation will set the given value as the mole fracti
 
 __Task B: Add output__:
 
-In the next step, we want to add some output to the simulation.
-First we want to know the water mass in the (porous-medium) system.
-Therefore, we evaluate the storage term of the water component.
+In the next step, we want to add some output to the simulation. The standard method for providing simulation output for visualization is via a `VtkWriter`. These tools take the grid geometries of each domain, as well as the solutions and write spatial output that one can view in visualization tools such as paraview.
+
+Although this Vtk output is very useful, some output is more suited for other forms of visualization.
+Two examples of data ouput formats are `.csv` files and `.json` files.
+From these files, data can be visualized using many different tools.
+In this exercise, data exported to a `.csv` file will be plotted using gnuplot, and data exported to a `.json` output will be plotted using the matplotlib python library.
+
+First as example output data, we want to investigate the water mass in the (porous-medium) system.
+There are two ways to evaluate this: (1) the total storage of water mass in the porous medium domain, and (2) the water vapor flux along the interface.
+The total storage can be plotted over time, whereas the water vapor flux will vary spatially along the interface as well as over time.
+
+First, we evaluate the storage term of the water component.
 * Have a look at the function `evaluateWaterMassStorageTerm()` in the `porousmediumsubproblem.hh`.
   Then implement a calculation of the total water mass:
   $`\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
-  Afterwards, adapt the method `init()` such that the variable `initialWaterContent_` is initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
-
-We also want to investigate the temporal evolution of the water mass.
-The following steps need to be done to do so.
-Check if all instructions are implemented accordingly:
-* Calculate the initial water mass at the beginning of the simulation and add the water mass loss to the output.
-  Based on the water mass loss you can derive the evaporation rate. Because despite at the interface, all
-  boundaries have Neumann no-flow conditions and no sources are present, the water can only leave the system
-  via the porous-medium free-flow interface. The evaporation in [mm/d] can be calculated by:
+  Afterwards, adapt the method `startStorageEvaluation()` such that the variable `initialWaterContent_` is initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
+
+In order to evaluate how the stored water mass evolves over time, we should make sure the following is implemented:
+* Calculate the initial water mass at the beginning of the simulation.
+* For each time step, add the water mass loss to the output.
+* Based on the water mass loss you can derive the evaporation rate. The evaporation in [mm/d] can be calculated by:
   $`e = \frac{\textrm{d}\, M^\text{w}}{\textrm{d}\, t} / A_\textrm{interface}`$.
+  (all boundaries have Neumann no-flow conditions and no sources are present, meaning the water can only leave the system
+  via the porous-medium free-flow interface.)
+
+During the simulation a `.csv` file is created (enable by setting `[Problem.ExportStorage] = true` in the `params.input` file). In addition, a gnuplot script `StorageOverTime.gp` is written and the mass loss and cumulative mass are plotted over time (`StorageOverTime.png`), when plotting is enabled with `[Problem.PlotStorage] = true`.
 
-Finally we want to know the distribution of the water mass fluxes across the interface.
+Second, we want to know the distribution of the water mass fluxes across the interface.
 * Have a look at the function `evaluateInterfaceFluxes()` in the porous medium problem.
   Use the facilities therein to return the values of ...massCouplingCondition... from the `couplingManager`
-  for each coupling scvf. Then the fluxes are visualized with gnuplot, when setting `Problem.PlotFluxes = true`.
-  If the simulation is too fast, you can have a look at the flux*.png files after the simulation.
-* You can use the parameter `Problem.PlotStorage = true` to see the temporal evolution of the evaporation rate
-  and the cumulative water mass loss.
+  for each coupling scvf. Multiply this with the relevant face sizes, extrusion factor, massfraction and seconds per day to get the [mm/d] units as seen in the storage evaluation.
+* When the `[Problem.ExportFluxes] = true` parameter is enabled, simulation data will be exported for this simulation to a `.json` file.
+  This file can flexibly handle data of different types, which in this case is helpful as we have both temporal and spatial data.
+* Use the python plotting script `plotFluxes.py` to visualize the flux distribution across the surface at different times.
+  This script uses `matplotlib`, a very popular python based visualization library.
 
 Compile and run the simulation and take a look at the results.
 
@@ -302,7 +327,7 @@ Now you are able to simulate a complete drying of the porous medium.
 
 Several RANS turbulence models are implemented in DuMu<sup>x</sup>.
 This part of the exercise consists of the following steps:
-* replacing the Navier-Stokes model by the zero equation turbulence model,
+* replacing the Navier-Stokes model by the K-Omega SST turbulence model,
 * switching to a symmetry boundary condition,
 * applying a grid refinement towards the interface,
 * subsequently refining the grid (convergence study).
@@ -317,46 +342,55 @@ The file `freeflowsubproblem.hh` is your free flow problem file  and `properties
 For using the compositional zero equation turbulence model, the following header files need to be included
  in properties file:
 ```
-#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/compositional/sstncmodel.hh>
 ```
 and in problem file:
 ```
-#include <dumux/freeflow/rans/problem.hh>
+#include <dumux/freeflow/turbulencemodel.hh>
+#include <dumux/freeflow/turbulenceproperties.hh>
+#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
 #include <dumux/freeflow/rans/boundarytypes.hh>
 ```
 
 The includes for the NavierStokesNC model and the NavierStokesProblem are no longer needed and can be removed.
 
 Make sure your free flow problem inherits from the correct parent type:
-* Change the entry in the `StokesZeroEq` definition accordingly (non-isothermal zero equation model, ZeroEqNCNI) in the properties file,
-* Adapt the inheritance of the problem class in problem file.
-
-Take a look into the two headers included above to see how the correct TypeTag and the Problem class the inherit from are called.
-
-Here, the turbulent free flow is wall bounded which means that the main reason for the development
-of turbulent flow is the presence of walls.
-The porous medium at the bottom and also the channel wall are such walls.
-Inside the problem you have to tell the turbulence model, where these walls are located
-by providing boundary conditions at the wall. For that, use
- ```cpp
+* Change the entry in the `FreeflowModel` definition accordingly (multi-component non-isothermal K-Omega SST model, SSTNCNI) in the properties file,
+* Adapt the inheritance of the problem class in problem file (Use `RANSProblem<TypeTag>` rather than `NavierStokesStaggeredProblem<TypeTag>`).
+  Take a look into the two headers included above to see how the correct TypeTag and the Problem class the inherit from are called.
+
+Here, the turbulent free-flow is wall-bounded, meaning shear stress and turbulence in the flow develop primarily due to the walls.
+The porous medium at the bottom and also the channel wall are examples of such walls.
+Within the problem the location of boundaries representing walls need to be defined.
+To do this, add the following function to the correct locations within the `boundaryTypes` function:
+```cpp
   values.setWall();
- ```
+```
 
-In addition, especially for the zero-equation models, any element in the free-flow domain interacts with the walls,
-e.g. this defines the wall distance which is needed to calculate the eddy viscosity.
-To get all these interactions, you have to call
- ```cpp
- stokesProblem->updateStaticWallProperties();
- ```
-in `main.cc`.
-However, there is also a solution-dependent component of these interactions, e.g. for a correct
-damping of the eddy viscosity toward the wall, the velocity gradient at the wall and inside the
-cells is needed.
-These dynamic interactions are to be updated by calling
+With the locations of the wall designated, and the problem initialized, a series of constant spatial properties,
+in particular the distance to the nearest wall from each cell center, should be initialized.
+To do this, add the following to the `main.cc` file.
 ```cpp
-stokesProblem->updateDynamicWallProperties(stokesSol);
+ freeflowProblem->updateStaticWallProperties();
 ```
-in the time loop (after `// Update dynamic wall properties`).
+
+In addition, there is also a non-local solution-dependent aspect of the turbulence models which is to be updated at the end of each step.
+An example of this is a stencil extended velocity gradient, and other flow properties in relation to the wall.
+These dynamic interactions are to be initialized in the mainfile directly after the `updateStaticWallProperties()` call,
+as well as in the time loop after `// Update dynamic wall properties`:
+```cpp
+freeflowProblem->updateDynamicWallProperties(freeflowSol);
+```
+
+In addition to designating the locations of walls,
+additional boundary conditions and initial conditions need to be set for the two new primary variables $k$ and $\omega$.
+In the `boundaryTypes` function, set both variables on all walls to be dirichlet, except for the right boundary, which should have outflow conditions.
+For the initial conditions, reynolds number specific base conditions should be applied everywhere.
+In the problem constructor, uncomment the code calcualting these terms,
+then apply the `turbulentKineticEnergy_`and `dissipation_` variables to their primary variables in all locations,
+except for on the wall boundaries, where these values can be set to zero.
+In addition, dirichlet constraints for the dissipation or $\omega$ variable will be set for all wall adjacent cells.
+This is done in the `isDirichletCell` function, as well as the `dirichlet` function already, and requires no further changes.
 
 Compile and run your new coupled problem and take a look at the results in Paraview.
 In addition to the standard variables and parameters, you can now analyze turbulence model specific quantities
@@ -368,13 +402,13 @@ The result for the turbulent viscosity should look like this:
 
 __Task B__:
 
-Instead of computing the whole cross-section of a channel, you can use symmetric boundary conditions at the top boundary of your free flow domain by replacing all previous boundary conditions with
+Instead of computing the whole cross-section of a channel,
+you can use symmetric boundary conditions at the top boundary of your free flow domain by replacing all previous boundary conditions at the top with
 ```c++
 values.setAllSymmetry();
 ```
 
-In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWallAtPos(globalPos)`
-and `initialAtPos(globalPos)` method.
+In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `initialAtPos(globalPos)` method.
 
 __Task C__:
 
@@ -393,11 +427,11 @@ __Task D__:
 
 For the grid convergence study, run various simulations with the following grading parameters:
 ```c++
-* [Stokes.Grid] Grading1 = 1.2, [Darcy.Grid] Grading1 = -1.2
-* [Stokes.Grid] Grading1 = 1.3, [Darcy.Grid] Grading1 = -1.3
-* [Stokes.Grid] Grading1 = 1.4, [Darcy.Grid] Grading1 = -1.4
-* [Stokes.Grid] Grading1 = 1.5, [Darcy.Grid] Grading1 = -1.5
-* [Stokes.Grid] Grading1 = 1.6, [Darcy.Grid] Grading1 = -1.6
+* [Freeflow.Grid] Grading1 = 1.2, [PorousMedium.Grid] Grading1 = -1.2
+* [Freeflow.Grid] Grading1 = 1.3, [PorousMedium.Grid] Grading1 = -1.3
+* [Freeflow.Grid] Grading1 = 1.4, [PorousMedium.Grid] Grading1 = -1.4
+* [Freeflow.Grid] Grading1 = 1.5, [PorousMedium.Grid] Grading1 = -1.5
+* [Freeflow.Grid] Grading1 = 1.6, [PorousMedium.Grid] Grading1 = -1.6
 ```
 
 By changing the parameter `Problem.Name` for each grading factor, you avoid losing the `.vtu` and `.pvd` files of the previous simulation runs.
diff --git a/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt
index be512539cde3ac912481c09798f81eef382e9eb0..748d97ec6644ce0abe5658149d2c1ab43a9e3d94 100644
--- a/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt
@@ -1,6 +1,7 @@
 # executables for ex_interface_coupling_ff-pm
 dumux_add_test(NAME exercise_interface_coupling_ff-pm
-               SOURCES main.cc)
+               SOURCES main.cc
+               LABELS ffpm)
 
 # add a symlink for each input file
 add_input_file_links()
diff --git a/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh b/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
index a7984532a58d3c52b3ed907e7e275ce32980d687..0f7d590fd107443694000c2b227642a5d187a0dc 100644
--- a/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief The free flow sub problem
  */
-#ifndef DUMUX_STOKES_SUBPROBLEM_HH
-#define DUMUX_STOKES_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_SUBPROBLEM_HH
 
 #include <dumux/freeflow/navierstokes/staggered/problem.hh>
 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
@@ -60,22 +60,15 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
     using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
 
 public:
-    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                       std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"),
+    eps_(1e-6),
+    couplingManager_(couplingManager)
     {
         deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
     }
 
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
-
-
     /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary segment.
@@ -163,17 +156,15 @@ public:
         return values;
     }
 
+    /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
-    //! Set the coupling manager
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
-
-    //! Get the coupling manager
-    const CouplingManager& couplingManager() const
-    { return *couplingManager_; }
-
-
-   /*!
+    /*!
      * \brief Evaluate the initial value for a control volume.
      *
      * \param globalPos The global position
@@ -192,17 +183,13 @@ public:
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().couplingData().darcyPermeability(element, scvf);
-    }
+    { return couplingManager().couplingData().darcyPermeability(element, scvf); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
     /*!
      * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
@@ -241,6 +228,14 @@ public:
         return analyticalVelocityX_;
     }
 
+    //! Set the coupling manager
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    //! Get the coupling manager
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
+
     // \}
 
 private:
diff --git a/exercises/exercise-coupling-ff-pm/interface/main.cc b/exercises/exercise-coupling-ff-pm/interface/main.cc
index 5cba8d6bd44c87c2a740dc31940b7001f94dc4fc..c7cce2fd18504e4452c88439b0c127a93f88489c 100644
--- a/exercises/exercise-coupling-ff-pm/interface/main.cc
+++ b/exercises/exercise-coupling-ff-pm/interface/main.cc
@@ -61,25 +61,25 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesOneP;
-    using DarcyTypeTag = Properties::TTag::DarcyOneP;
+    using FreeflowTypeTag = Properties::TTag::FreeflowOneP;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumFlowOneP;
 
     //TODO: dumux-course-task 1.C
     // ******************** comment-out this section for the last exercise **************** //
 
     // create two individual grids (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 
     // ************************************************************************************ //
 
@@ -104,13 +104,13 @@ int main(int argc, char** argv)
 //
 //     Params params;
 //
-//     auto elementSelectorStokes = [&](const auto& element)
+//     auto elementSelectorFreeflow = [&](const auto& element)
 //     {
 //         double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
 //         return element.geometry().center()[1] > interface;
 //     };
 //
-//     auto elementSelectorDarcy = [&](const auto& element)
+//     auto elementSelectorPorousMedium = [&](const auto& element)
 //     {
 //         double interface  =  params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
 //         return element.geometry().center()[1] < interface;
@@ -118,92 +118,96 @@ int main(int argc, char** argv)
 //
 //     using SubGrid = Dune::SubGrid<dim, HostGrid>;
 //
-//     Dumux::GridManager<SubGrid> subGridManagerStokes;
-//     Dumux::GridManager<SubGrid> subGridManagerDarcy;
+//     Dumux::GridManager<SubGrid> subGridManagerFreeflow;
+//     Dumux::GridManager<SubGrid> subGridManagerPorousMedium;
 //
 //     // initialize subgrids
-//     subGridManagerStokes.init(hostGrid, elementSelectorStokes, "Stokes");
-//     subGridManagerDarcy.init(hostGrid, elementSelectorDarcy, "Darcy");
+//     subGridManagerFreeflow.init(hostGrid, elementSelectorFreeflow, "Freeflow");
+//     subGridManagerPorousMedium.init(hostGrid, elementSelectorPorousMedium, "PorousMedium");
 //
 //     // we compute on the leaf grid view
-//     const auto& darcyGridView = subGridManagerDarcy.grid().leafGridView();
-//     const auto& stokesGridView = subGridManagerStokes.grid().leafGridView();
+//     const auto& porousMediumGridView = subGridManagerPorousMedium.grid().leafGridView();
+//     const auto& freeflowGridView = subGridManagerFreeflow.grid().leafGridView();
 
     // ************************************************************************************ //
 
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
 
-    // initialize the vtk output module
-    const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+    // intialize the vtk output module
+    const auto freeflowName = getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
 
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
 
     //TODO: dumux-course-task 1.B
     //****** uncomment the add analytical solution of v_x *****//
-    // stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+    // freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
+    freeflowVtkWriter.write(0.0);
 
-    stokesVtkWriter.write(0.0);
+    using PorousMediumSolutionVector = GetPropType<PorousMediumTypeTag, Properties::SolutionVector>;
+    VtkOutputModule<PorousMediumGridVariables, PorousMediumSolutionVector> porousMediumVtkWriter(*porousMediumGridVariables,
+                                                                                                 sol[porousMediumIdx],
+                                                                                                 porousMediumName);
 
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
 
     // the assembler for a stationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager);
 
     // the linear solver
@@ -218,8 +222,8 @@ int main(int argc, char** argv)
     nonLinearSolver.solve(sol);
 
     // write vtk output
-    stokesVtkWriter.write(1.0);
-    darcyVtkWriter.write(1.0);
+    freeflowVtkWriter.write(1.0);
+    porousMediumVtkWriter.write(1.0);
 
     Parameters::print();
 
diff --git a/exercises/exercise-coupling-ff-pm/interface/params.input b/exercises/exercise-coupling-ff-pm/interface/params.input
index 22730c83706e84a4c14db95473405ca48c0be3bd..eab4c116545ff15c89fd02e2ea78671d5d31b49a 100644
--- a/exercises/exercise-coupling-ff-pm/interface/params.input
+++ b/exercises/exercise-coupling-ff-pm/interface/params.input
@@ -1,5 +1,5 @@
- # for dune-subgrid
- [Grid]
+# for dune-subgrid
+[Grid]
 Positions0 = 0 1
 Positions1 = 0 0.2 0.3 0.65
 Cells0 = 100
@@ -9,7 +9,7 @@ Amplitude = 0.04 # [m]
 Offset = 0.5 # [m]
 Scaling = 0.2 #[m]
 
-[Stokes.Grid]
+[Freeflow.Grid]
 Verbosity = true
 Positions0 = 0.0 1.0
 Positions1 = 1.0 2.0
@@ -17,7 +17,7 @@ Cells0 = 20
 Cells1 = 100
 Grading1 = 1
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 Verbosity = true
 Positions0 = 0.0 1.0
 Positions1 = 0.0 1.0
@@ -25,13 +25,13 @@ Cells0 = 20
 Cells1 = 20
 Grading1 = 1
 
-[Stokes.Problem]
-Name = stokes
+[Freeflow.Problem]
+Name = freeflow
 PressureDifference = 1e-9
 EnableInertiaTerms = false
 
-[Darcy.Problem]
-Name = darcy
+[PorousMedium.Problem]
+Name = porousmedium
 
 [SpatialParams]
 Permeability = 1e-6 # m^2
diff --git a/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh b/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
index 836e041cbb39f40d35fd89b4deeaf6effb3170a3..68967766e1998018df981fbf5236772c7943d702 100644
--- a/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
@@ -21,8 +21,8 @@
  *
  * \brief The porous medium flow sub problem
  */
-#ifndef DUMUX_DARCY_SUBPROBLEM_HH
-#define DUMUX_DARCY_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
 
 #include <dumux/porousmediumflow/problem.hh>
 #include <dumux/common/properties.hh>
@@ -60,18 +60,19 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
-                   std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+                           std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"),
+    eps_(1e-7),
+    couplingManager_(couplingManager)
     {}
 
-
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
     {
         BoundaryTypes values;
@@ -90,7 +91,7 @@ public:
         return values;
     }
 
-        /*!
+    /*!
      * \brief Evaluate the boundary conditions for a Dirichlet control volume.
      *
      * \param element The element for which the Dirichlet boundary condition is set
@@ -150,8 +151,6 @@ public:
                        const SubControlVolume &scv) const
     { return NumEqVector(0.0); }
 
-    // \}
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -161,11 +160,7 @@ public:
      * variables.
      */
     PrimaryVariables initial(const Element &element) const
-    {
-        return PrimaryVariables(0.0);
-    }
-
-    // \}
+    { return PrimaryVariables(0.0); }
 
     //! Set the coupling manager
     void setCouplingManager(std::shared_ptr<CouplingManager> cm)
diff --git a/exercises/exercise-coupling-ff-pm/interface/properties.hh b/exercises/exercise-coupling-ff-pm/interface/properties.hh
index b1b82bbd158b8f77f84f25e01321dfdda4b9b02b..3f5d388d1f8df46f952774fe9fb4e54a19f9fef0 100644
--- a/exercises/exercise-coupling-ff-pm/interface/properties.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/properties.hh
@@ -52,39 +52,39 @@ namespace Dumux::Properties {
 
 // Create new type tags
 namespace TTag {
-struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
-struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct FreeflowOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct PorousMediumFlowOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesOneP>
+struct CouplingManager<TypeTag, TTag::FreeflowOneP>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumFlowOneP>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyOneP>
+struct CouplingManager<TypeTag, TTag::PorousMediumFlowOneP>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowOneP, Properties::TTag::FreeflowOneP, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumFlowOneP> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowOneP> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 
 // the fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyOneP>
+struct FluidSystem<TypeTag, TTag::PorousMediumFlowOneP>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
 };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesOneP>
+struct FluidSystem<TypeTag, TTag::FreeflowOneP>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
@@ -92,7 +92,7 @@ struct FluidSystem<TypeTag, TTag::StokesOneP>
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyOneP>
+struct Grid<TypeTag, TTag::PorousMediumFlowOneP>
 {
     static constexpr auto dim = 2;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -107,7 +107,7 @@ struct Grid<TypeTag, TTag::DarcyOneP>
     // using type = Dune::SubGrid<dim, HostGrid>;
 };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesOneP>
+struct Grid<TypeTag, TTag::FreeflowOneP>
 {
     static constexpr auto dim = 2;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -123,16 +123,16 @@ struct Grid<TypeTag, TTag::StokesOneP>
 };
 
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyOneP> {
+struct SpatialParams<TypeTag, TTag::PorousMediumFlowOneP> {
     using type = OnePSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
 };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 
 } //end namespace Dumux::Properties
 
diff --git a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
index 8edaf9f1b6abc8c44db23b26a695e1b16a81bb27..b44a8c7a6e6c961c81bfcc9d0949ede28ce31d4e 100644
--- a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -1,6 +1,6 @@
 # executables for ex_interface_coupling_ff-pm
 dumux_add_test(NAME exercise_models_coupling_ff-pm
-               SOURCES main.cc)
+               SOURCES main.cc
+               LABELS ffpm)
 
-# add a symlink for each input file
-add_input_file_links()
+dune_symlink_to_source_files(FILES "params.input" plotFluxes.py)
diff --git a/exercises/exercise-coupling-ff-pm/models/freeflowsubproblem.hh b/exercises/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
index 57d749f2b05a1be9e324c66b18c7cf4b9181ead3..227e276af325321fe70ddf435172519bf84893c0 100644
--- a/exercises/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
@@ -21,8 +21,8 @@
  * \ingroup NavierStokesTests
  * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
  */
-#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
-#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -69,24 +69,17 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
     static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
-    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                       std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"),
+    eps_(1e-6),
+    couplingManager_(couplingManager)
     {
         velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
         pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
         moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
     }
 
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
-
-
     /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary segment.
@@ -173,21 +166,7 @@ public:
         return values;
     }
 
-
     /*!
-     * \brief Set the coupling manager
-     */
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
-
-    /*!
-     * \brief Get the coupling manager
-     */
-    const CouplingManager& couplingManager() const
-    { return *couplingManager_; }
-
-
-   /*!
      * \brief Evaluate the initial value for a control volume.
      *
      * \param globalPos The global position
@@ -216,26 +195,40 @@ public:
         return values;
     }
 
-    void setTimeLoop(TimeLoopPtr timeLoop)
-    { timeLoop_ = timeLoop; }
+    /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().couplingData().darcyPermeability(element, scvf);
-    }
+    { return couplingManager().couplingData().darcyPermeability(element, scvf); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
-    // \}
+    /*!
+     * \brief Set the coupling manager
+     */
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    /*!
+     * \brief Get the coupling manager
+     */
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
+
+    void setTimeLoop(TimeLoopPtr timeLoop)
+    { timeLoop_ = timeLoop; }
 
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -271,6 +264,7 @@ private:
         const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
         fluidState.setEnthalpy(0, enthalpy);
     }
+
     // the height of the free-flow domain
     const Scalar height_() const
     { return this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1]; }
diff --git a/exercises/exercise-coupling-ff-pm/models/main.cc b/exercises/exercise-coupling-ff-pm/models/main.cc
index 2fe38c8e2497454b3e6d5661d3687d79fd992d2e..736d7bec9027711719ab5821dc3fecc8b226c0f1 100644
--- a/exercises/exercise-coupling-ff-pm/models/main.cc
+++ b/exercises/exercise-coupling-ff-pm/models/main.cc
@@ -19,7 +19,7 @@
 /*!
  * \file
  *
- * \brief A test problem for the coupled Stokes/Darcy problem (1p)
+ * \brief A test problem for the coupled Freeflow/PorousMedium problem (1p)
  */
 #include <config.h>
 
@@ -61,51 +61,51 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesNC;
-    using DarcyTypeTag = Properties::TTag::DarcyOnePNC;
+    using FreeflowTypeTag = Properties::TTag::FreeflowNC;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumOnePNC;
 
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+    GetPropType<FreeflowTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+    using Scalar = GetPropType<FreeflowTypeTag, Properties::Scalar>;
     const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
     const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -113,58 +113,61 @@ int main(int argc, char** argv)
     // instantiate time loop
     auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
-    stokesProblem->setTimeLoop(timeLoop);
-    darcyProblem->setTimeLoop(timeLoop);
+    freeflowProblem->setTimeLoop(timeLoop);
+    porousMediumProblem->setTimeLoop(timeLoop);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
     auto solOld = sol;
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
-
-    // initialize the vtk output module
-    const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
-
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
-    stokesVtkWriter.write(0.0);
-
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
-
-    // initialize the subproblems
-    darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
+
+    // intialize the vtk output module
+    const auto freeflowName = getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
+
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
+    freeflowVtkWriter.write(0.0);
+
+    VtkOutputModule<PorousMediumGridVariables, GetPropType<PorousMediumTypeTag, Properties::SolutionVector>> porousMediumVtkWriter(*porousMediumGridVariables,
+                                                                                                                                   sol[porousMediumIdx], porousMediumName);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
+
+    // Start the water mass storage calculation
+    porousMediumProblem->startStorageEvaluation(sol[porousMediumIdx], *porousMediumGridVariables);
+    // Start the evaporation flux calculation
+    porousMediumProblem->startFluxEvaluation();
 
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager,
                                                  timeLoop, solOld);
 
@@ -186,20 +189,20 @@ int main(int argc, char** argv)
 
         // make the new solution the old solution
         solOld = sol;
-        stokesGridVariables->advanceTimeStep();
-        darcyGridVariables->advanceTimeStep();
+        freeflowGridVariables->advanceTimeStep();
+        porousMediumGridVariables->advanceTimeStep();
 
         // advance to the time loop to the next step
         timeLoop->advanceTimeStep();
 
         // call the postTimeStep routine for output
-        darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+        porousMediumProblem->postTimeStep(sol[porousMediumIdx], *porousMediumGridVariables);
 
         // write vtk output
         if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
         {
-            stokesVtkWriter.write(timeLoop->time());
-            darcyVtkWriter.write(timeLoop->time());
+            freeflowVtkWriter.write(timeLoop->time());
+            porousMediumVtkWriter.write(timeLoop->time());
         }
 
         // report statistics of this time step
@@ -210,8 +213,8 @@ int main(int argc, char** argv)
 
     } while (!timeLoop->finished());
 
-    timeLoop->finalize(stokesGridView.comm());
-    timeLoop->finalize(darcyGridView.comm());
+    timeLoop->finalize(freeflowGridView.comm());
+    timeLoop->finalize(porousMediumGridView.comm());
 
     Parameters::print();
 
diff --git a/exercises/exercise-coupling-ff-pm/models/params.input b/exercises/exercise-coupling-ff-pm/models/params.input
index 152a908fcd026862e0353843f42fbf456cc8d016..1fa6c1c23a65f5233cd451ba9b8acf5d822a6b15 100644
--- a/exercises/exercise-coupling-ff-pm/models/params.input
+++ b/exercises/exercise-coupling-ff-pm/models/params.input
@@ -3,25 +3,25 @@ DtInitial = 1 # s
 EpisodeLength = -360 # s # 0.25 days
 TEnd = 256000 # s # 2 days
 
-[Stokes.Grid]
+[Freeflow.Grid]
 LowerLeft = 0 1
 UpperRight = 1 2
 Cells = 16 16
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 UpperRight = 1 1
 Cells = 16 16
 
-[Stokes.Problem]
-Name = stokes
+[Freeflow.Problem]
+Name = freeflow
 EnableGravity = true
 EnableInertiaTerms = true
 MoleFraction = 0.0 # -
 Pressure = 1e5 # Pa
 Velocity = 1 # m/s
 
-[Darcy.Problem]
-Name = darcy
+[PorousMedium.Problem]
+Name = porousmedium
 EnableGravity = true
 Saturation = 0.1 # -
 MoleFraction = 0.1 # -
@@ -42,8 +42,9 @@ Temperature = 293.15
 
 [Problem]
 Name = models_coupling
-PlotFluxes = false
-PlotStorage = false
+ExportStorage = true
+PlotStorage = true
+ExportFluxes = true
 
 [Newton]
 MaxSteps = 12
diff --git a/exercises/exercise-coupling-ff-pm/models/plotFluxes.py b/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
new file mode 100644
index 0000000000000000000000000000000000000000..832062f0fc2b4c9945f22d1c384cb413a2f7d151
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
@@ -0,0 +1,32 @@
+import json
+import matplotlib.pyplot as plt
+
+# JSON file
+jsonFilePath = "./flux_models_coupling.json"
+file = open(jsonFilePath)
+fullEvaporationData = json.loads(file.read())
+
+for key in fullEvaporationData:
+    print("developing plot " + key)
+
+    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
+    plt.title('Interfacial Evaporation Rate: ' + key, fontsize=24)
+    plt.xlabel('Interface x')
+    plt.ylabel('Evaporation flux')
+    evaporationData = fullEvaporationData[key]
+    coordinates = evaporationData["FaceCenterCoordinates"]
+    timeSeries = evaporationData["Time"]
+    xCoords = [x for x,y in coordinates]
+    facewiseEvaporationResults = evaporationData["FaceEvaporation"]
+
+    count = 0
+    for time, results in zip(timeSeries, facewiseEvaporationResults):
+        if (count%2 == 0 and time >= 1.0): # plot every 5th time step
+            ax.plot(xCoords, results, label="{:0.3f}".format(time))
+        count = count+1
+
+    ax.legend(title="Time (s)")
+    outputFig = plt.gcf()
+    outputFig.savefig('./EvaporationOverInterface_' + key + '.png', bbox_inches='tight')
+    plt.clf()
+    plt.close()
diff --git a/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh b/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
index 01d64fdd8ca7064b080958e04a126dee25c0c5db..b05d07338875d6084c8d70c257a3cc7866096719 100644
--- a/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
@@ -21,14 +21,15 @@
  *
  * \brief A simple Darcy test problem (cell-centered finite volume method).
  */
-#ifndef DUMUX_DARCY_SUBPROBLEM_HH
-#define DUMUX_DARCY_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
 #include <dumux/common/timeloop.hh>
 #include <dumux/common/numeqvector.hh>
 #include <dumux/io/gnuplotinterface.hh>
+#include <dumux/common/metadata.hh>
 
 #include <dumux/porousmediumflow/problem.hh>
 
@@ -76,15 +77,34 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
-                   std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+                           std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"),
+    eps_(1e-7),
+    couplingManager_(couplingManager)
     {
         moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
 
-        // initialize output file
-        plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+        // initialize storage output file (.csv)
+        exportStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportStorage", false);
         plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
-        storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+        if (exportStorage_)
+        {
+            storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+            initializeStorageOutput();
+        }
+
+        // initialize flux output file (.json)
+        exportFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportFluxes", false);
+        if (exportFluxes_)
+        {
+            fluxFileName_ = "flux" + getParam<std::string>("Problem.Name");
+            simulationKey_ = getParam<std::string>("Problem.Name", "case1");
+            initializeFluxOutput();
+        }
+    }
+
+    void initializeStorageOutput()
+    {
         storageFile_.open(storageFileName_);
         storageFile_ << "#Time[s]" << ";"
                      << "WaterMass[kg]" << ";"
@@ -93,19 +113,56 @@ public:
                      << std::endl;
     }
 
+    void initializeFluxOutput()
+    {
+        Dumux::MetaData::Collector collector;
+        if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+            Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+        Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+    }
 
-    /*!
-     * \brief Initialize the problem.
-     */
     template<class SolutionVector, class GridVariables>
-    void init(const SolutionVector& curSol,
-              const GridVariables& gridVariables)
+    void startStorageEvaluation(const SolutionVector& curSol,
+                                const GridVariables& gridVariables)
     {
         // TODO: dumux-course-task 2.B
         // Initialize `initialWaterContent_` and assign that to `lastWaterMass_`.
 
     }
 
+    void startFluxEvaluation()
+    {
+        // Get Interface Data
+        std::vector<GlobalPosition> faceLocation;
+        std::vector<Scalar> faceArea;
+
+        for (const auto& element : elements(this->gridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->gridGeometry());
+            fvGeometry.bindElement(element);
+
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+                    continue;
+                faceLocation.push_back(scvf.ipGlobal());
+                faceArea.push_back(scvf.area());
+            }
+        }
+        std::vector<Scalar> faceEvaporation(faceArea.size(), 0.0);
+
+        Dumux::MetaData::Collector collector;
+        if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+            Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+        collector[simulationKey_]["Time"]= {0.0};
+        collector[simulationKey_]["TimeStepIdx"]= {0};
+        collector[simulationKey_]["FaceCenterCoordinates"] = faceLocation;
+        collector[simulationKey_]["FaceLength"] = faceArea;
+        collector[simulationKey_]["FaceEvaporation"] = {faceEvaporation};
+        collector[simulationKey_]["FluxOverFullSurface"] = {0.0};
+        Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+    }
+
     template<class SolutionVector, class GridVariables>
     void postTimeStep(const SolutionVector& curSol,
                       const GridVariables& gridVariables)
@@ -114,7 +171,7 @@ public:
         evaluateWaterMassStorageTerm(curSol, gridVariables);
         evaluateInterfaceFluxes(curSol, gridVariables);
 
-        if (plotStorage_)
+        if (exportStorage_ && plotStorage_)
         {
             gnuplotStorage_.resetPlot();
             gnuplotStorage_.setDatafileSeparator(';');
@@ -127,7 +184,7 @@ public:
             gnuplotStorage_.setOption("set y2range [0.0:0.5]");
             gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
             gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
-            gnuplotStorage_.plot("temp");
+            gnuplotStorage_.plot("StorageOverTime");
         }
     }
 
@@ -163,11 +220,14 @@ public:
                                  / timeLoop_->timeStepSize();
         lastWaterMass_ = waterMass;
 
-        storageFile_ << timeLoop_->time() << ";"
-                     << waterMass << ";"
-                     << cumMassLoss << ";"
-                     << evaporationRate
-                     << std::endl;
+        if (exportStorage_)
+        {
+            storageFile_ << timeLoop_->time() << ";"
+                         << waterMass << ";"
+                         << cumMassLoss << ";"
+                         << evaporationRate
+                         << std::endl;
+        }
 
         return waterMass;
     }
@@ -177,8 +237,7 @@ public:
                                  const GridVariables& gridVariables)
 
     {
-        std::vector<Scalar> x;
-        std::vector<Scalar> y;
+        std::vector<Scalar> faceEvaporation;
 
         for (const auto& element : elements(this->gridGeometry().gridView()))
         {
@@ -196,35 +255,32 @@ public:
                 // TODO: dumux-course-task 2.B
                 // Use "massCouplingCondition" from the couplingManager here
                 NumEqVector flux(0.0);
-
-                x.push_back(scvf.center()[0]);
-                y.push_back(flux[transportCompIdx]);
+                faceEvaporation.push_back(flux[transportCompIdx]);
             }
         }
+        Scalar fluxOverFullSurface = std::accumulate(faceEvaporation.begin(), faceEvaporation.end(), 0.0);
 
-        if (plotFluxes_)
+        if (exportFluxes_)
         {
-            gnuplotInterfaceFluxes_.resetPlot();
-            gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
-            gnuplotInterfaceFluxes_.setXRange(this->gridGeometry().bBoxMin()[0], this->gridGeometry().bBoxMax()[0]);
-            gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
-            gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
-            gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
-            std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
-                                       "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
-            gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
-            gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+            Dumux::MetaData::Collector collector;
+            if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+                Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+            collector[simulationKey_]["Time"].push_back(time());
+            collector[simulationKey_]["TimeStepIdx"].push_back(timeStepIndex());
+            collector[simulationKey_]["FaceEvaporation"].push_back(faceEvaporation);
+            collector[simulationKey_]["FluxOverFullSurface"].push_back(fluxOverFullSurface);
+            Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+            std::cout << "Total flux per local method = " << fluxOverFullSurface << "\n";
         }
     }
 
-
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
     {
         BoundaryTypes values;
@@ -277,7 +333,6 @@ public:
                        const SubControlVolume& scv) const
     { return NumEqVector(0.0); }
 
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -288,20 +343,18 @@ public:
      */
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
-        static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+        static const Scalar freeflowPressure = getParamFromGroup<Scalar>("Freeflow", "Problem.Pressure");
 
         PrimaryVariables values(0.0);
 
         // TODO: dumux-course-task 2.A
         // Declare here which phases are present.
 
-        values[Indices::pressureIdx] = stokesPressure;
+        values[Indices::pressureIdx] = freeflowPressure;
         values[transportCompIdx] = moleFraction_;
         return values;
     }
 
-    // \}
-
     //! Set the coupling manager
     void setCouplingManager(std::shared_ptr<CouplingManager> cm)
     { couplingManager_ = cm; }
@@ -313,6 +366,12 @@ public:
     void setTimeLoop(TimeLoopPtr timeLoop)
     { timeLoop_ = timeLoop; }
 
+    Scalar time() const
+    { return timeLoop_->time(); }
+
+    int timeStepIndex() const
+    { return timeLoop_->timeStepIndex(); }
+
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
     { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
@@ -328,19 +387,22 @@ private:
 
     Scalar eps_;
     Scalar moleFraction_;
-
-    Scalar initialWaterContent_ = 0.0;
-    Scalar lastWaterMass_ = 0.0;
-
     TimeLoopPtr timeLoop_;
     std::shared_ptr<CouplingManager> couplingManager_;
 
+    // Storage calculation and output
+    Scalar initialWaterContent_ = 0.0;
+    Scalar lastWaterMass_ = 0.0;
     std::string storageFileName_;
     std::ofstream storageFile_;
-    bool plotFluxes_;
+    bool exportStorage_;
     bool plotStorage_;
-    Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
     Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+
+    // Flux output
+    std::string fluxFileName_;
+    std::string simulationKey_;
+    bool exportFluxes_;
 };
 
 } //end namespace Dumux
diff --git a/exercises/exercise-coupling-ff-pm/models/properties.hh b/exercises/exercise-coupling-ff-pm/models/properties.hh
index 72b63dabf59ac4ae386451daa05a04e68f8617fc..d3313c732a9bfbcf0211f5e9f8582710c8a59061 100644
--- a/exercises/exercise-coupling-ff-pm/models/properties.hh
+++ b/exercises/exercise-coupling-ff-pm/models/properties.hh
@@ -57,39 +57,39 @@ namespace Dumux::Properties {
 namespace TTag {
 // TODO: dumux-course-task 2.A
 // Change to property of the `FluidSystem` such that `H2OAir` is used directly.
-struct DarcyOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
-struct StokesNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+struct PorousMediumOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+struct FreeflowNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesNC>
+struct CouplingManager<TypeTag, TTag::FreeflowNC>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePNC>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumOnePNC>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyOnePNC>
+struct CouplingManager<TypeTag, TTag::PorousMediumOnePNC>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesNC, Properties::TTag::StokesNC, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowNC, Properties::TTag::FreeflowNC, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyOnePNC> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumOnePNC> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesNC> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowNC> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 
 // The fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyOnePNC>
+struct FluidSystem<TypeTag, TTag::PorousMediumOnePNC>
 {
     using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesNC>
+struct FluidSystem<TypeTag, TTag::FreeflowNC>
 {
     using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
@@ -97,25 +97,25 @@ struct FluidSystem<TypeTag, TTag::StokesNC>
 
 // Use moles
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::DarcyOnePNC> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::PorousMediumOnePNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 
 // Do not replace one equation with a total mass balance
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePNC> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::PorousMediumOnePNC> { static constexpr int value = 3; };
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::StokesNC> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::FreeflowNC> { static constexpr int value = 3; };
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyOnePNC> { using type = Dune::YaspGrid<2>; };
+struct Grid<TypeTag, TTag::PorousMediumOnePNC> { using type = Dune::YaspGrid<2>; };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::FreeflowNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 //! Use a model with constant tortuosity for the effective diffusivity
 template<class TypeTag>
-struct EffectiveDiffusivityModel<TypeTag, TTag::DarcyOnePNC>
+struct EffectiveDiffusivityModel<TypeTag, TTag::PorousMediumOnePNC>
 { using type = DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>; };
 
 // TODO: dumux-course-task 2.A
@@ -125,16 +125,16 @@ struct EffectiveDiffusivityModel<TypeTag, TTag::DarcyOnePNC>
 template<class TypeTag>
 // TODO: dumux-course-task 2.A
 // Adapt the spatial params here.
-struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+struct SpatialParams<TypeTag, TTag::PorousMediumOnePNC> {
     using type = OnePSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
 };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 
 } // end namespace Dumux::Properties
 
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
index b36e7939761d666cf2edd664c1094c524791acaf..d669f666b71c240af706a173c1fa5de01aa191a0 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -1,6 +1,7 @@
 # executables for ex_interface_coupling_ff-pm
 dumux_add_test(NAME exercise_turbulence_coupling_ff-pm
-               SOURCES main.cc)
+               SOURCES main.cc
+               LABELS ffpm)
 
 # add a symlink for each input file
 add_input_file_links()
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
index 4887b0962d62c58ef2b32dddc4245dceca01ad7b..5a68f51303acfc7ef060cce07447f35855f5f69a 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief The free-flow sub problem
  */
-#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
-#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -46,6 +46,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
     using ParentType = NavierStokesStaggeredProblem<TypeTag>;
 
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+    static constexpr auto dimWorld = GridView::dimensionworld;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -56,6 +57,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
 
     using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
@@ -75,8 +77,11 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
     static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
-    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                       std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"),
+    eps_(1e-6),
+    couplingManager_(couplingManager)
     {
         refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
         refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
@@ -84,19 +89,25 @@ public:
         refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
 
         diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
-                                                                     getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+                                                                     getParamFromGroup<std::string>(this->paramGroup(),
+                                                                     "Problem.InterfaceDiffusionCoefficientAvg"));
+
+// TODO:   // ******************** uncomment this section for the first task 3.A ****************** //
+//         FluidSystem::init();
+//         Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
+//         FluidState fluidState;
+//         const auto phaseIdx = 0;
+//         fluidState.setPressure(phaseIdx, refPressure_);
+//         fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
+//         fluidState.setMassFraction(phaseIdx, phaseIdx, 1.0);
+//         Scalar density = FluidSystem::density(fluidState, phaseIdx);
+//         Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, phaseIdx) / density;
+//         Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+//         turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(refVelocity_, diameter, kinematicViscosity);
+//         dissipation_ = turbulenceProperties.dissipationRate(refVelocity_, diameter, kinematicViscosity);
+//         // ************************************************************************************* //
     }
 
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
-
-
     /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary segment.
@@ -120,7 +131,8 @@ public:
         }
 
 // TODO: dumux-course-task 3.A
-// set boundary conditions for the turbulence model at the wall (values.setWall()) at the upper and lower boundary
+// set wall conditions for the turbulence model at the wall (values.setWall()) at the upper and lower boundary
+// set boundary conditions for the turbulence model primary variables everywhere (outflow on right boundary, otherwise dirichlet)
         if (onLowerBoundary_(globalPos))
         {
             values.setDirichlet(Indices::velocityXIdx);
@@ -160,16 +172,65 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+     * \brief Evaluate the boundary conditions for dirichlet values at the boundary.
      *
-     * \param element The element
-     * \param scvf The subcontrolvolume face
+     * \param element The finite element
+     * \param scvf the sub control volume face
+     * \note used for cell-centered discretization schemes
      */
-    PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+    PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
-        values = initialAtPos(pos);
+        const auto globalPos = scvf.ipGlobal();
+        PrimaryVariables values(initialAtPos(globalPos));
+        return values;
+    }
 
+    /*!
+     * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param scv The sub control volume
+     * \param pvIdx The primary variable index in the solution vector
+     */
+    template<class Element, class FVElementGeometry, class SubControlVolume>
+    bool isDirichletCell(const Element& element,
+                         const FVElementGeometry& fvGeometry,
+                         const SubControlVolume& scv,
+                         int pvIdx) const
+    {
+        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+        {
+            // For the komega model we set a fixed dissipation (omega) for all cells at the wall
+            for (const auto& scvf : scvfs(fvGeometry))
+                if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx)
+                    return true;
+        }
+        return false;
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for fixed values at cell centers
+     *
+     * \param element The finite element
+     * \param scv the sub control volume
+     * \note used for cell-centered discretization schemes
+     */
+    PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
+    {
+        const auto globalPos = scv.center();
+        PrimaryVariables values(initialAtPos(globalPos));
+        unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
+
+        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+        {
+            // For the komega model we set a fixed value for the dissipation
+            const auto wallDistance = ParentType::wallDistance(elementIdx);
+            using std::pow;
+            values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
+                                                    / (ParentType::betaOmega() * wallDistance * wallDistance);
+            return values;
+        }
         return values;
     }
 
@@ -201,21 +262,15 @@ public:
         return values;
     }
 
-
     /*!
-     * \brief Set the coupling manager
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
      */
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
-     * \brief Get the coupling manager
-     */
-    const CouplingManager& couplingManager() const
-    { return *couplingManager_; }
-
-
-   /*!
      * \brief Evaluate the initial value for a control volume.
      *
      * \param globalPos The global position
@@ -233,6 +288,18 @@ public:
         values[Indices::velocityXIdx] = refVelocity();
         values[Indices::temperatureIdx] = refTemperature();
 
+        // TODO: dumux-course-task 3.A
+        // Set initial conditions for the TKE and the Dissipation. Values calculated in the constructor
+//         values[Indices::turbulentKineticEnergyIdx] = TODO??;
+//         values[Indices::dissipationIdx] = TODO??;
+//         // TODO: dumux-course-task 3.B
+//         // Remove the condition `onUpperBoundary_(globalPos)` here.
+//         if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+//         {
+//             values[Indices::turbulentKineticEnergyIdx] = 0.0;
+//             values[Indices::dissipationIdx] = 0.0;
+//         }
+
         // TODO: dumux-course-task 3.B
         // Remove the condition `onUpperBoundary_(globalPos)` here.
         if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
@@ -257,7 +324,6 @@ public:
     const Scalar refTemperature() const
     { return refTemperature_; }
 
-
     void setTimeLoop(TimeLoopPtr timeLoop)
     { timeLoop_ = timeLoop; }
 
@@ -265,19 +331,25 @@ public:
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center()); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
-    // \}
+    /*!
+     * \brief Set the coupling manager
+     */
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    /*!
+     * \brief Get the coupling manager
+     */
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
 
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -325,6 +397,8 @@ private:
     Scalar refPressure_;
     Scalar refMoleFrac_;
     Scalar refTemperature_;
+    Scalar turbulentKineticEnergy_;
+    Scalar dissipation_;
 
     TimeLoopPtr timeLoop_;
 
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/main.cc b/exercises/exercise-coupling-ff-pm/turbulence/main.cc
index 769e1f4f78af35ce582cbbfbac1f3d8f60c659ff..cfe3ffeda883992b0f49d21643e2ba5ba57ea842 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/main.cc
+++ b/exercises/exercise-coupling-ff-pm/turbulence/main.cc
@@ -62,51 +62,51 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesZeroEq;
-    using DarcyTypeTag = Properties::TTag::DarcyTwoPTwoCNI;
+    using FreeflowTypeTag = Properties::TTag::FreeflowModel;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumFlowModel;
 
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+    GetPropType<FreeflowTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+    using Scalar = GetPropType<FreeflowTypeTag, Properties::Scalar>;
     const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
     const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -116,23 +116,23 @@ int main(int argc, char** argv)
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // set timeloop for the subproblems, needed for boundary value variations
-    stokesProblem->setTimeLoop(timeLoop);
-    darcyProblem->setTimeLoop(timeLoop);
+    freeflowProblem->setTimeLoop(timeLoop);
+    porousMediumProblem->setTimeLoop(timeLoop);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
     auto solOld = sol;
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
 
     // TODO: dumux-course-task 3.A
     // Update static wall properties
@@ -141,36 +141,37 @@ int main(int argc, char** argv)
     // Update dynamic wall properties
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
-
-    // initialize the vtk output module
-    const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
-
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
-    stokesVtkWriter.write(0.0);
-
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
+
+    // intialize the vtk output module
+    const auto freeflowName = getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
+
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
+    freeflowVtkWriter.write(0.0);
+
+    using PorousMediumSolutionVector = GetPropType<PorousMediumTypeTag, Properties::SolutionVector>;
+    VtkOutputModule<PorousMediumGridVariables, PorousMediumSolutionVector> porousMediumVtkWriter(*porousMediumGridVariables, sol[porousMediumIdx], porousMediumName);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
 
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager,
                                                  timeLoop, solOld);
 
@@ -194,19 +195,19 @@ int main(int argc, char** argv)
         // TODO: dumux-course-task 3.A
         // Update dynamic wall properties
 
-        // post time step treatment of Darcy problem
-        darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+        // post time step treatment of PorousMedium problem
+        porousMediumProblem->postTimeStep(sol[porousMediumIdx], *porousMediumGridVariables, timeLoop->timeStepSize());
 
         // advance grid variables to the next time step
-        stokesGridVariables->advanceTimeStep();
-        darcyGridVariables->advanceTimeStep();
+        freeflowGridVariables->advanceTimeStep();
+        porousMediumGridVariables->advanceTimeStep();
 
         // advance to the time loop to the next step
         timeLoop->advanceTimeStep();
 
         // write vtk output
-        stokesVtkWriter.write(timeLoop->time());
-        darcyVtkWriter.write(timeLoop->time());
+        freeflowVtkWriter.write(timeLoop->time());
+        porousMediumVtkWriter.write(timeLoop->time());
 
         // report statistics of this time step
         timeLoop->reportTimeStep();
@@ -216,8 +217,8 @@ int main(int argc, char** argv)
 
     } while (!timeLoop->finished());
 
-    timeLoop->finalize(stokesGridView.comm());
-    timeLoop->finalize(darcyGridView.comm());
+    timeLoop->finalize(freeflowGridView.comm());
+    timeLoop->finalize(porousMediumGridView.comm());
 
     Parameters::print();
 
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/params.input b/exercises/exercise-coupling-ff-pm/turbulence/params.input
index 7f0d58a81b9fbb723c235c7d3c085034f23b3154..b64b055899ecfe2729bc27c74946883b02d496b5 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/params.input
+++ b/exercises/exercise-coupling-ff-pm/turbulence/params.input
@@ -3,7 +3,7 @@ DtInitial =  1e-1 # [s]
 MaxTimeStepSize = 43200 # [s] (12 hours)
 TEnd = 864000 # [s] (6 days)
 
-[Stokes.Grid]
+[Freeflow.Grid]
 Positions0 = 0.0 0.25
 Positions1 = 0.25 0.5
 Grading0 = 1.0
@@ -12,7 +12,7 @@ Cells0 = 15
 Cells1 = 20
 Verbosity = true
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 Positions0 = 0.0 0.25
 Positions1 = 0.0 0.25
 Cells0 = 15
@@ -21,16 +21,16 @@ Grading0 = 1.0
 Grading1 = 1.0
 Verbosity = true
 
-[Stokes.Problem]
-Name = stokes
+[Freeflow.Problem]
+Name = freeflow
 RefVelocity = 3.5 # [m/s]
 RefPressure = 1e5 # [Pa]
 refMoleFrac = 0 # [-]
 RefTemperature = 298.15 # [K]
 EnableInertiaTerms = true
 
-[Darcy.Problem]
-Name = darcy
+[PorousMedium.Problem]
+Name = porousmedium
 Pressure = 1e5
 Saturation = 0.5 # initial Sw
 Temperature = 298.15 # [K]
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
index 5d989c367192c6fe81776a07768e540357fe08c3..bc2bd75a63131f1ac416822e15b1ce2e5858ec7a 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
@@ -21,8 +21,8 @@
  *
  * \brief The porous medium sub problem
  */
-#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
-#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW2P2C_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW2P2C_SUBPROBLEM_HH
 
 #include <dumux/porousmediumflow/problem.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -77,8 +77,10 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
-                   std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+                           std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"),
+    eps_(1e-7),
+    couplingManager_(couplingManager)
     {
         pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
         initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
@@ -89,7 +91,6 @@ public:
                                                                      getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
     }
 
-
     template<class SolutionVector, class GridVariables>
     void postTimeStep(const SolutionVector& curSol,
                       const GridVariables& gridVariables,
@@ -123,12 +124,12 @@ public:
     }
 
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
     {
         BoundaryTypes values;
@@ -207,8 +208,6 @@ public:
                        const SubControlVolume &scv) const
     { return NumEqVector(0.0); }
 
-    // \}
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -227,8 +226,6 @@ public:
         return values;
     }
 
-    // \}
-
     /*!
      * \brief Set the coupling manager
      */
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/properties.hh b/exercises/exercise-coupling-ff-pm/turbulence/properties.hh
index 3d86bc5d12bf8207f61ffad6c292f9ae6e9f1d85..14bc29ff7ff023b927eb5ad26b808e6d8fbf04be 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/properties.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/properties.hh
@@ -51,43 +51,43 @@ namespace Dumux::Properties {
 
 // Create new type tags
 namespace TTag {
-struct DarcyTwoPTwoCNI { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+struct PorousMediumFlowModel { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
 // TODO: dumux-course-task 3.A
-// Change the entry in the `StokesZeroEq` definition accordingly.
-struct StokesZeroEq { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+// Change the entry in the `FreeflowModel` definition accordingly.
+struct FreeflowModel { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesZeroEq>
+struct CouplingManager<TypeTag, TTag::FreeflowModel>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCNI>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumFlowModel>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct CouplingManager<TypeTag, TTag::PorousMediumFlowModel>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesZeroEq, Properties::TTag::StokesZeroEq, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowModel, Properties::TTag::FreeflowModel, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumFlowModel> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesZeroEq> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowModel> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::PorousMediumFlowModel> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesZeroEq> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::FreeflowModel> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // the fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
+struct FluidSystem<TypeTag, TTag::PorousMediumFlowModel> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesZeroEq>
+struct FluidSystem<TypeTag, TTag::FreeflowModel>
 {
   using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
   static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
@@ -96,32 +96,32 @@ struct FluidSystem<TypeTag, TTag::StokesZeroEq>
 
 // Use formulation based on mole fractions
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::PorousMediumFlowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 
 // The gas component balance (air) is replaced by the total mass balance
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::PorousMediumFlowModel> { static constexpr int value = 3; };
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::StokesZeroEq> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::FreeflowModel> { static constexpr int value = 3; };
 
 
 //! Set the default formulation to pw-Sn: This can be over written in the problem.
 template<class TypeTag>
-struct Formulation<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct Formulation<TypeTag, TTag::PorousMediumFlowModel>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct SpatialParams<TypeTag, TTag::PorousMediumFlowModel>
 { using type = TwoPSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>; };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 
 } // end namespace Dumux::Properties
 
diff --git a/exercises/extradoc/ex_ff-pm-turb_diffusivity.png b/exercises/extradoc/ex_ff-pm-turb_diffusivity.png
index c7173d0e3bfc0146d77f5a4401005cf9c7df215b..796c223d467c9ab722c2128dae7976f5ea4c1015 100644
Binary files a/exercises/extradoc/ex_ff-pm-turb_diffusivity.png and b/exercises/extradoc/ex_ff-pm-turb_diffusivity.png differ
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
index 1a5c23ded8e44237fc8ae9dbc01888c14e04e57d..8ef01b7e68d413ea597e4ac2956ca6647abf8835 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
@@ -1,18 +1,27 @@
 dumux_add_test(NAME exercise_interface_coupling_ff-pm_original
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=0)
 
 dumux_add_test(NAME exercise_interface_coupling_ff-pm_a_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=1)
 
 dumux_add_test(NAME exercise_interface_coupling_ff-pm_b_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=2)
 
 dumux_add_test(NAME exercise_interface_coupling_ff-pm_c_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=3)
 
+dumux_add_test(NAME exercise_interface_coupling_ff-pm_extra_solution
+               SOURCES main.cc
+               LABELS ffpm
+               COMPILE_DEFINITIONS EXNUMBER=4)
+
 # add a symlink for each input file
 add_input_file_links()
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
index 2438cadae30a813c3e5c6fa65a360f832bfd6a5d..e083d38aece298468b187e93ee673462738703ae 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
@@ -20,8 +20,8 @@
  * \file
  * \brief The free flow sub problem
  */
-#ifndef DUMUX_STOKES_SUBPROBLEM_HH
-#define DUMUX_STOKES_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_INTERFACE_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_INTERFACE_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -62,21 +62,11 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
 
 public:
     FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"), eps_(1e-6), couplingManager_(couplingManager)
     {
         deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
     }
 
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
-
-
     /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary segment.
@@ -91,7 +81,7 @@ public:
 
         const auto& globalPos = scvf.dofPosition();
 
-#if EXNUMBER == 0 // flow from top to bottom
+#if EXNUMBER == 0  // flow from top to bottom
         if(onUpperBoundary_(globalPos))
         {
             values.setDirichlet(Indices::velocityXIdx);
@@ -103,7 +93,7 @@ public:
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
         }
-#else // flow flom left to right
+#elif EXNUMBER < 4 // flow flom left to right
         if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
             values.setDirichlet(Indices::pressureIdx);
         else
@@ -111,6 +101,14 @@ public:
             values.setDirichlet(Indices::velocityXIdx);
             values.setDirichlet(Indices::velocityYIdx);
         }
+#else
+        if (onRightBoundary_(globalPos))
+            values.setDirichlet(Indices::pressureIdx);
+        else
+        {
+            values.setDirichlet(Indices::velocityXIdx);
+            values.setDirichlet(Indices::velocityYIdx);
+        }
 #endif
 
         if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
@@ -191,8 +189,15 @@ public:
     const CouplingManager& couplingManager() const
     { return *couplingManager_; }
 
+    /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
-   /*!
+    /*!
      * \brief Evaluate the initial value for a control volume.
      *
      * \param globalPos The global position
@@ -202,6 +207,9 @@ public:
         PrimaryVariables values(0.0);
 #if EXNUMBER == 0
         values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
+#elif EXNUMBER == 4
+        values[Indices::velocityXIdx] = 1e-6 * (globalPos[1] - this->gridGeometry().bBoxMin()[1])
+                                      * (this->gridGeometry().bBoxMax()[1] - globalPos[1]);
 #else
         // set fixed pressures on the left and right boundary
         if(onLeftBoundary_(globalPos))
@@ -209,7 +217,6 @@ public:
         if(onRightBoundary_(globalPos))
             values[Indices::pressureIdx] = 0.0;
 #endif
-
         return values;
     }
 
@@ -217,17 +224,13 @@ public:
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().couplingData().darcyPermeability(element, scvf);
-    }
+    { return couplingManager().couplingData().darcyPermeability(element, scvf); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
     /*!
      * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/main.cc b/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
index 73426d92eae6dc2f7db1e3066dc1a96e77aa7aca..b4db7eb11730c8da40b60fc4f839535cde3baefe 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
@@ -60,23 +60,23 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesOneP;
-    using DarcyTypeTag = Properties::TTag::DarcyOneP;
+    using FreeflowTypeTag = Properties::TTag::FreeflowOneP;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumOneP;
 
 #if EXNUMBER < 3
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 #else
     // use dune-subgrid to create the individual grids
     static constexpr int dim = 2;
@@ -96,13 +96,13 @@ int main(int argc, char** argv)
 
     Params params;
 
-    auto elementSelectorStokes = [&](const auto& element)
+    auto elementSelectorFreeflow = [&](const auto& element)
     {
         double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
         return element.geometry().center()[1] > interface;
     };
 
-    auto elementSelectorDarcy = [&](const auto& element)
+    auto elementSelectorPorousMedium = [&](const auto& element)
     {
         double interface  =  params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
         return element.geometry().center()[1] < interface;
@@ -110,91 +110,92 @@ int main(int argc, char** argv)
 
     using SubGrid = Dune::SubGrid<dim, HostGrid>;
 
-    Dumux::GridManager<SubGrid> subGridManagerStokes;
-    Dumux::GridManager<SubGrid> subGridManagerDarcy;
+    Dumux::GridManager<SubGrid> subGridManagerFreeflow;
+    Dumux::GridManager<SubGrid> subGridManagerPorousMedium;
 
     // initialize subgrids
-    subGridManagerStokes.init(hostGrid, elementSelectorStokes, "Stokes");
-    subGridManagerDarcy.init(hostGrid, elementSelectorDarcy, "Darcy");
+    subGridManagerFreeflow.init(hostGrid, elementSelectorFreeflow, "Freeflow");
+    subGridManagerPorousMedium.init(hostGrid, elementSelectorPorousMedium, "PorousMedium");
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = subGridManagerDarcy.grid().leafGridView();
-    const auto& stokesGridView = subGridManagerStokes.grid().leafGridView();
+    const auto& porousMediumGridView = subGridManagerPorousMedium.grid().leafGridView();
+    const auto& freeflowGridView = subGridManagerFreeflow.grid().leafGridView();
 
 #endif
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
 
-    // initialize the vtk output module
-    const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
 
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+    // intialize the vtk output module
+    const auto freeflowName = getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
+
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
 
 #if EXNUMBER >= 2
-    stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+    freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
 #endif
 
-    stokesVtkWriter.write(0.0);
+    freeflowVtkWriter.write(0.0);
 
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
+    VtkOutputModule<PorousMediumGridVariables, GetPropType<PorousMediumTypeTag, Properties::SolutionVector>> porousMediumVtkWriter(*porousMediumGridVariables, sol[porousMediumIdx], porousMediumName);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
 
     // the assembler for a stationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager);
 
     // the linear solver
@@ -209,8 +210,8 @@ int main(int argc, char** argv)
     nonLinearSolver.solve(sol);
 
     // write vtk output
-    stokesVtkWriter.write(1.0);
-    darcyVtkWriter.write(1.0);
+    freeflowVtkWriter.write(1.0);
+    porousMediumVtkWriter.write(1.0);
 
     Parameters::print();
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/params.input b/exercises/solution/exercise-coupling-ff-pm/interface/params.input
index 6935e9946b93abd3c903b08b32c3cb15151c6b84..2f1dc08b438a052ce1ef5ff5ab57b34612c2c991 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/params.input
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/params.input
@@ -1,5 +1,5 @@
- # for dune-subgrid
- [Grid]
+# for dune-subgrid
+[Grid]
 Positions0 = 0 1
 Positions1 = 0 0.2 0.3 0.65
 Cells0 = 100
@@ -9,7 +9,7 @@ Amplitude = 0.04 # [m]
 Offset = 0.5 # [m]
 Scaling = 0.2 #[m]
 
-[Stokes.Grid]
+[Freeflow.Grid]
 Verbosity = true
 Positions0 = 0.0 1.0
 Positions1 = 1.0 2.0
@@ -17,7 +17,7 @@ Cells0 = 20
 Cells1 = 100
 Grading1 = 1
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 Verbosity = true
 Positions0 = 0.0 1.0
 Positions1 = 0.0 1.0
@@ -25,12 +25,12 @@ Cells0 = 20
 Cells1 = 20
 Grading1 = 1
 
-[Stokes.Problem]
+[Freeflow.Problem]
 Name = stokes
 PressureDifference = 1e-9
 EnableInertiaTerms = false
 
-[Darcy.Problem]
+[PorousMedium.Problem]
 Name = darcy
 
 [SpatialParams]
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
index 173c5c9d5a0b20f5fb8ce097d75ff6b445206434..85e1978c2d7fe20a8785d431c2d2cfce8ddafc31 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
@@ -21,8 +21,8 @@
  *
  * \brief The porous medium flow sub problem
  */
-#ifndef DUMUX_DARCY_SUBPROBLEM_HH
-#define DUMUX_DARCY_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW_INTERFACE_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW_INTERFACE_SUBPROBLEM_HH
 
 #include <dumux/porousmediumflow/problem.hh>
 #include <dumux/common/properties.hh>
@@ -61,17 +61,16 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
                    std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"), eps_(1e-7), couplingManager_(couplingManager)
     {}
 
-
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
     {
         BoundaryTypes values;
@@ -88,7 +87,7 @@ public:
         return values;
     }
 
-        /*!
+    /*!
      * \brief Evaluate the boundary conditions for a Dirichlet control volume.
      *
      * \param element The element for which the Dirichlet boundary condition is set
@@ -145,8 +144,6 @@ public:
                        const SubControlVolume &scv) const
     { return NumEqVector(0.0); }
 
-    // \}
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -156,11 +153,7 @@ public:
      * variables.
      */
     PrimaryVariables initial(const Element &element) const
-    {
-        return PrimaryVariables(0.0);
-    }
-
-    // \}
+    { return PrimaryVariables(0.0);}
 
     //! Set the coupling manager
     void setCouplingManager(std::shared_ptr<CouplingManager> cm)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh b/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
index 642c815f4cc4c9261b12f7430e7cb68e6af046dd..f8a8587c2f9c1299ae570b97d53ed1755142a8e1 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
@@ -52,33 +52,33 @@ namespace Dumux::Properties {
 
 // Create new type tags
 namespace TTag {
-struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
-struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct FreeflowOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+struct PorousMediumOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesOneP>
+struct CouplingManager<TypeTag, TTag::FreeflowOneP>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumOneP>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyOneP>
+struct CouplingManager<TypeTag, TTag::PorousMediumOneP>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowOneP, Properties::TTag::FreeflowOneP, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumOneP> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowOneP> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyOneP>
+struct Grid<TypeTag, TTag::PorousMediumOneP>
 {
     static constexpr auto dim = 2;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -92,7 +92,7 @@ struct Grid<TypeTag, TTag::DarcyOneP>
 #endif
 };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesOneP>
+struct Grid<TypeTag, TTag::FreeflowOneP>
 {
     static constexpr auto dim = 2;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -108,29 +108,29 @@ struct Grid<TypeTag, TTag::StokesOneP>
 
 // the fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyOneP>
+struct FluidSystem<TypeTag, TTag::PorousMediumOneP>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
 };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesOneP>
+struct FluidSystem<TypeTag, TTag::FreeflowOneP>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
 };
 
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyOneP> {
+struct SpatialParams<TypeTag, TTag::PorousMediumOneP> {
     using type = OnePSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
 };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
 
 } // end namespace Dumux::Properties
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
index 7b21c78003a37e90903339b58870ca0a1da47100..a7f4ba40ce22825867d7b9c7a10cc4f5548053b5 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -1,26 +1,29 @@
 dumux_add_test(NAME exercise_models_coupling_ff-pm_original
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=0
                COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_original
-               CMD_ARGS -Problem.PlotStorage 0 -Problem.PlotFluxes 0)
+               CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
 
 dumux_add_test(NAME exercise_models_coupling_ff-pm_a_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=1
                COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_a_solution
-               CMD_ARGS -Problem.PlotStorage 0 -Problem.PlotFluxes 0)
+               CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
 
 dumux_add_test(NAME exercise_models_coupling_ff-pm_b_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=2
                COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_b_solution
-               CMD_ARGS -Problem.PlotStorage 0 -Problem.PlotFluxes 0)
+               CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
 
 dumux_add_test(NAME exercise_models_coupling_ff-pm_c_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=3
                COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_c_solution
-               CMD_ARGS -Problem.PlotStorage 0 -Problem.PlotFluxes 0)
+               CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
 
-# add a symlink for each input file
-add_input_file_links()
+dune_symlink_to_source_files(FILES "params.input" plotFluxes.py)
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/freeflowsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
index fca504852d8783857b8b50bc90fadee323ed6c45..d8ae90cde7d62a9915c673a355dc7731e7301f31 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/freeflowsubproblem.hh
@@ -21,8 +21,8 @@
  * \ingroup NavierStokesTests
  * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
  */
-#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
-#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_MODEL_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_MODEL_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -73,23 +73,13 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
 
 public:
     FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"), eps_(1e-6), couplingManager_(couplingManager)
     {
         velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
         pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
         moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
     }
 
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
-
-
     /*!
      * \brief Specifies which kind of boundary condition should be
      *        used for which equation on a given boundary segment.
@@ -176,25 +166,11 @@ public:
         return values;
     }
 
-
-    /*!
-     * \brief Set the coupling manager
-     */
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
-
     /*!
-     * \brief Get the coupling manager
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param globalPos The global position
      */
-    const CouplingManager& couplingManager() const
-    { return *couplingManager_; }
-
-
-    /*!
-      * \brief Evaluate the initial value for a control volume.
-      *
-      * \param globalPos The global position
-      */
      PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
      {
          PrimaryVariables values(0.0);
@@ -219,6 +195,14 @@ public:
          return values;
      }
 
+    /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
+
     void setTimeLoop(TimeLoopPtr timeLoop)
     { timeLoop_ = timeLoop; }
 
@@ -226,19 +210,25 @@ public:
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().couplingData().darcyPermeability(element, scvf);
-    }
+    { return couplingManager().couplingData().darcyPermeability(element, scvf); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
-    // \}
+    /*!
+     * \brief Set the coupling manager
+     */
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    /*!
+     * \brief Get the coupling manager
+     */
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
 
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/main.cc b/exercises/solution/exercise-coupling-ff-pm/models/main.cc
index f1d4b67a1c57f766d413af2343f42e3dddceb703..6416f2b7e9dc7432d1d0f4926a92f1d791892fec 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/main.cc
+++ b/exercises/solution/exercise-coupling-ff-pm/models/main.cc
@@ -54,9 +54,6 @@ int main(int argc, char** argv)
 {
     using namespace Dumux;
 
-    // initialize MPI, finalize is done automatically on exit
-    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
-
     // initialize MPI+x, finalize is done automatically on exit
     Dumux::initialize(argc, argv);
 
@@ -64,51 +61,51 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesNC;
-    using DarcyTypeTag = Properties::TTag::DarcyOnePNC;
+    using FreeflowTypeTag = Properties::TTag::FreeflowNC;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumOnePNC;
 
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+    GetPropType<FreeflowTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+    using Scalar = GetPropType<FreeflowTypeTag, Properties::Scalar>;
     const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
     const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -116,64 +113,66 @@ int main(int argc, char** argv)
     // instantiate time loop
     auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
-    stokesProblem->setTimeLoop(timeLoop);
-    darcyProblem->setTimeLoop(timeLoop);
+    freeflowProblem->setTimeLoop(timeLoop);
+    porousMediumProblem->setTimeLoop(timeLoop);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
     auto solOld = sol;
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
 
     // initialize the vtk output module
 #if EXNUMBER >= 1
     const std::array<std::string, 3> part = {"a", "b", "c"};
-    const auto stokesName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+    const auto freeflowName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
 #else
-    const auto stokesName = "orig_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = "orig_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+    const auto freeflowName = "orig_" + getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = "orig_" + getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
 #endif
 
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
-    stokesVtkWriter.write(0);
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
+    freeflowVtkWriter.write(0.0);
 
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
+    VtkOutputModule<PorousMediumGridVariables, GetPropType<PorousMediumTypeTag, Properties::SolutionVector>> porousMediumVtkWriter(*porousMediumGridVariables, sol[porousMediumIdx], porousMediumName);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
 
-    // initialize the subproblems
-    darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+    // Start the water mass storage calculation
+    porousMediumProblem->startStorageEvaluation(sol[porousMediumIdx], *porousMediumGridVariables);
+    // Start the evaporation flux calculation
+    porousMediumProblem->startFluxEvaluation();
 
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager,
                                                  timeLoop, solOld);
 
@@ -195,20 +194,20 @@ int main(int argc, char** argv)
 
         // make the new solution the old solution
         solOld = sol;
-        stokesGridVariables->advanceTimeStep();
-        darcyGridVariables->advanceTimeStep();
+        freeflowGridVariables->advanceTimeStep();
+        porousMediumGridVariables->advanceTimeStep();
 
         // advance to the time loop to the next step
         timeLoop->advanceTimeStep();
 
         // call the postTimeStep routine for output
-        darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+        porousMediumProblem->postTimeStep(sol[porousMediumIdx], *porousMediumGridVariables);
 
         // write vtk output
         if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
         {
-            stokesVtkWriter.write(timeLoop->time());
-            darcyVtkWriter.write(timeLoop->time());
+            freeflowVtkWriter.write(timeLoop->time());
+            porousMediumVtkWriter.write(timeLoop->time());
         }
 
         // report statistics of this time step
@@ -219,8 +218,8 @@ int main(int argc, char** argv)
 
     } while (!timeLoop->finished());
 
-    timeLoop->finalize(stokesGridView.comm());
-    timeLoop->finalize(darcyGridView.comm());
+    timeLoop->finalize(freeflowGridView.comm());
+    timeLoop->finalize(porousMediumGridView.comm());
 
     Parameters::print();
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/params.input b/exercises/solution/exercise-coupling-ff-pm/models/params.input
index d2c30bd98bbd476126938af001f4b4c733e19a41..dc2a61319d61d8739eb10201b4446330dff26a86 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/params.input
+++ b/exercises/solution/exercise-coupling-ff-pm/models/params.input
@@ -3,25 +3,25 @@ DtInitial = 100 # s
 EpisodeLength = -360 # s # 0.25 days
 TEnd = 256000 # s # 2 days
 
-[Stokes.Grid]
+[Freeflow.Grid]
 LowerLeft = 0 1
 UpperRight = 1 2
 Cells = 16 16
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 UpperRight = 1 1
 Cells = 16 16
 
-[Stokes.Problem]
-Name = stokes
+[Freeflow.Problem]
+Name = freeflow
 EnableGravity = true
 EnableInertiaTerms = true
 MoleFraction = 0.0 # -
 Pressure = 1e5 # Pa
 Velocity = 1 # m/s
 
-[Darcy.Problem]
-Name = darcy
+[PorousMedium.Problem]
+Name = porousmedium
 EnableGravity = true
 Saturation = 0.1 # -
 MoleFraction = 0.1 # -
@@ -42,9 +42,10 @@ PorousMediumTemperature = 293.15
 Temperature = 293.15
 
 [Problem]
-Name = ex_models_coupling
-PlotFluxes = true
+Name = models_coupling
+ExportStorage = true
 PlotStorage = true
+ExportFluxes = true
 
 [Newton]
 MaxSteps = 12
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
new file mode 120000
index 0000000000000000000000000000000000000000..667b97bf502ed5041690750bc3b7c16923ee6dd8
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
@@ -0,0 +1 @@
+../../../exercise-coupling-ff-pm/models/plotFluxes.py
\ No newline at end of file
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
index e49f63bbc13d6e9863b31863fb9b42a5f5390314..5f1b26ab2661ac2888af422eb5aa004e693af0a2 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
@@ -21,15 +21,16 @@
  *
  * \brief A simple Darcy test problem (cell-centered finite volume method).
  */
-#ifndef DUMUX_DARCY_SUBPROBLEM_HH
-#define DUMUX_DARCY_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
 #include <dumux/common/timeloop.hh>
 #include <dumux/common/numeqvector.hh>
-
 #include <dumux/io/gnuplotinterface.hh>
+#include <dumux/common/metadata.hh>
+
 #include <dumux/porousmediumflow/problem.hh>
 
 namespace Dumux {
@@ -81,8 +82,10 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
-                   std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+                           std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"),
+    eps_(1e-7),
+    couplingManager_(couplingManager)
     {
 #if EXNUMBER >= 3
         saturation_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
@@ -90,10 +93,27 @@ public:
         moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
 #endif
 
-        // initialize output file
-        plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+        // initialize storage output file (.csv)
+        exportStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportStorage", false);
         plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
-        storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+        if (exportStorage_)
+        {
+            storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+            initializeStorageOutput();
+        }
+
+        // initialize flux output file (.json)
+        exportFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportFluxes", false);
+        if (exportFluxes_)
+        {
+            fluxFileName_ = "flux_" + getParam<std::string>("Problem.Name");
+            simulationKey_ = getParam<std::string>("Problem.Name", "case1");
+            initializeFluxOutput();
+        }
+    }
+
+    void initializeStorageOutput()
+    {
         storageFile_.open(storageFileName_);
         storageFile_ << "#Time[s]" << ";"
                      << "WaterMass[kg]" << ";"
@@ -102,13 +122,17 @@ public:
                      << std::endl;
     }
 
+    void initializeFluxOutput()
+    {
+        Dumux::MetaData::Collector collector;
+        if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+            Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+        Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+    }
 
-    /*!
-     * \brief Initialize the problem.
-     */
     template<class SolutionVector, class GridVariables>
-    void init(const SolutionVector& curSol,
-              const GridVariables& gridVariables)
+    void startStorageEvaluation(const SolutionVector& curSol,
+                                const GridVariables& gridVariables)
     {
 #if EXNUMBER >= 2
         initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables);
@@ -116,6 +140,39 @@ public:
 #endif
     }
 
+    void startFluxEvaluation()
+    {
+        // Get Interface Data
+        std::vector<GlobalPosition> faceLocation;
+        std::vector<Scalar> faceArea;
+
+        for (const auto& element : elements(this->gridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->gridGeometry());
+            fvGeometry.bindElement(element);
+
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+                    continue;
+                faceLocation.push_back(scvf.ipGlobal());
+                faceArea.push_back(scvf.area());
+            }
+        }
+        std::vector<Scalar> faceEvaporation(faceArea.size(), 0.0);
+
+        Dumux::MetaData::Collector collector;
+        if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+            Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+        collector[simulationKey_]["Time"]= {0.0};
+        collector[simulationKey_]["TimeStepIdx"]= {0};
+        collector[simulationKey_]["FaceCenterCoordinates"] = faceLocation;
+        collector[simulationKey_]["FaceLength"] = faceArea;
+        collector[simulationKey_]["FaceEvaporation"] = {faceEvaporation};
+        collector[simulationKey_]["FluxOverFullSurface"] = {0.0};
+        Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+    }
+
     template<class SolutionVector, class GridVariables>
     void postTimeStep(const SolutionVector& curSol,
                       const GridVariables& gridVariables)
@@ -124,7 +181,7 @@ public:
         evaluateWaterMassStorageTerm(curSol, gridVariables);
         evaluateInterfaceFluxes(curSol, gridVariables);
 
-        if (plotFluxes_)
+        if (exportStorage_ && plotStorage_)
         {
             gnuplotStorage_.resetPlot();
             gnuplotStorage_.setDatafileSeparator(';');
@@ -137,7 +194,7 @@ public:
             gnuplotStorage_.setOption("set y2range [0.0:0.5]");
             gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
             gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
-            gnuplotStorage_.plot("temp");
+            gnuplotStorage_.plot("StorageOverTime");
         }
     }
 
@@ -183,11 +240,14 @@ public:
                                  / timeLoop_->timeStepSize();
         lastWaterMass_ = waterMass;
 
-        storageFile_ << timeLoop_->time() << ";"
-                     << waterMass << ";"
-                     << cumMassLoss << ";"
-                     << evaporationRate
-                     << std::endl;
+        if (exportStorage_)
+        {
+            storageFile_ << timeLoop_->time() << ";"
+                         << waterMass << ";"
+                         << cumMassLoss << ";"
+                         << evaporationRate
+                         << std::endl;
+        }
 
         return waterMass;
     }
@@ -197,8 +257,7 @@ public:
                                  const GridVariables& gridVariables)
 
     {
-        std::vector<Scalar> x;
-        std::vector<Scalar> y;
+        std::vector<Scalar> faceEvaporation;
 
         for (const auto& element : elements(this->gridGeometry().gridView()))
         {
@@ -214,39 +273,37 @@ public:
                     continue;
 
 #if EXNUMBER >= 2
-                NumEqVector flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
+                NumEqVector flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf)
+                                   * scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor() * FluidSystem::molarMass(1) * -1.0 * 86400.0;
 #else
                 NumEqVector flux(0.0); // add "massCouplingCondition" from the couplingManager here
 #endif
-
-                x.push_back(scvf.center()[0]);
-                y.push_back(flux[transportCompIdx]);
+                faceEvaporation.push_back(flux[transportCompIdx]);
             }
         }
+        Scalar fluxOverFullSurface = std::accumulate(faceEvaporation.begin(), faceEvaporation.end(), 0.0);
 
-        if (plotFluxes_)
+        if (exportFluxes_)
         {
-            gnuplotInterfaceFluxes_.resetPlot();
-            gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
-            gnuplotInterfaceFluxes_.setXRange(this->gridGeometry().bBoxMin()[0], this->gridGeometry().bBoxMax()[0]);
-            gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
-            gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
-            gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
-            std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
-                                       "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
-            gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
-            gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+            Dumux::MetaData::Collector collector;
+            if (Dumux::MetaData::jsonFileExists(fluxFileName_))
+                Dumux::MetaData::readJsonFile(collector, fluxFileName_);
+            collector[simulationKey_]["Time"].push_back(time());
+            collector[simulationKey_]["TimeStepIdx"].push_back(timeStepIndex());
+            collector[simulationKey_]["FaceEvaporation"].push_back(faceEvaporation);
+            collector[simulationKey_]["FluxOverFullSurface"].push_back(fluxOverFullSurface);
+            Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
+            std::cout << "Total flux per local method = " << fluxOverFullSurface << "\n";
         }
     }
 
-
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
     {
         BoundaryTypes values;
@@ -299,7 +356,6 @@ public:
                        const SubControlVolume& scv) const
     { return NumEqVector(0.0); }
 
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -310,7 +366,7 @@ public:
      */
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
-        static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+        static const Scalar freeflowPressure = getParamFromGroup<Scalar>("Freeflow", "Problem.Pressure");
 
         PrimaryVariables values(0.0);
 #if EXNUMBER >= 3
@@ -322,12 +378,10 @@ public:
 #else
         values[transportCompIdx] = moleFraction_;
 #endif
-        values[pressureIdx] = stokesPressure;
+        values[pressureIdx] = freeflowPressure;
         return values;
     }
 
-    // \}
-
     //! Set the coupling manager
     void setCouplingManager(std::shared_ptr<CouplingManager> cm)
     { couplingManager_ = cm; }
@@ -339,6 +393,12 @@ public:
     void setTimeLoop(TimeLoopPtr timeLoop)
     { timeLoop_ = timeLoop; }
 
+    Scalar time() const
+    { return timeLoop_->time(); }
+
+    int timeStepIndex() const
+    { return timeLoop_->timeStepIndex(); }
+
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
     { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
@@ -359,18 +419,22 @@ private:
     Scalar moleFraction_;
 #endif
 
-    Scalar initialWaterContent_ = 0.0;
-    Scalar lastWaterMass_ = 0.0;
-
     TimeLoopPtr timeLoop_;
     std::shared_ptr<CouplingManager> couplingManager_;
 
+    // Storage calculation and output
+    Scalar initialWaterContent_ = 0.0;
+    Scalar lastWaterMass_ = 0.0;
     std::string storageFileName_;
     std::ofstream storageFile_;
-    bool plotFluxes_;
+    bool exportStorage_;
     bool plotStorage_;
-    Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
     Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+
+    // Flux output
+    std::string fluxFileName_;
+    std::string simulationKey_;
+    bool exportFluxes_;
 };
 
 } //end namespace Dumux
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/properties.hh b/exercises/solution/exercise-coupling-ff-pm/models/properties.hh
index aa66bfb1565a86c24f7d33b6326fa4d4942cef40..ca80958014b76f08cce119f27791ab30af5d2838 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/properties.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/properties.hh
@@ -57,49 +57,49 @@ namespace Dumux::Properties {
 
 // Create new type tags
 namespace TTag {
-struct StokesNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+struct FreeflowNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
 #if EXNUMBER >= 1
-struct DarcyOnePNC { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
+struct PorousMediumOnePNC { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
 #else
-struct DarcyOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+struct PorousMediumOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
 #endif
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesNC>
+struct CouplingManager<TypeTag, TTag::FreeflowNC>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePNC>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumOnePNC>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyOnePNC>
+struct CouplingManager<TypeTag, TTag::PorousMediumOnePNC>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesNC, Properties::TTag::StokesNC, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowNC, Properties::TTag::FreeflowNC, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesNC> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowNC> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyOnePNC> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumOnePNC> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::FreeflowNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyOnePNC> { using type = Dune::YaspGrid<2>; };
+struct Grid<TypeTag, TTag::PorousMediumOnePNC> { using type = Dune::YaspGrid<2>; };
 
 // The fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesNC>
+struct FluidSystem<TypeTag, TTag::FreeflowNC>
 {
     using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
     using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
 };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyOnePNC>
+struct FluidSystem<TypeTag, TTag::PorousMediumOnePNC>
 {
     using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
 #if EXNUMBER == 0
@@ -111,47 +111,47 @@ struct FluidSystem<TypeTag, TTag::DarcyOnePNC>
 
 // Do not replace one equation with a total mass balance
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::StokesNC> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::FreeflowNC> { static constexpr int value = 3; };
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePNC> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::PorousMediumOnePNC> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::DarcyOnePNC> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::PorousMediumOnePNC> { static constexpr bool value = true; };
 
 #if EXNUMBER >= 1
 //! Set the default formulation to pw-Sn: This can be over written in the problem.
 template<class TypeTag>
-struct Formulation<TypeTag, TTag::DarcyOnePNC>
+struct Formulation<TypeTag, TTag::PorousMediumOnePNC>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 #endif
 
 // Set the spatial parameters type
 #if EXNUMBER >= 1
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+struct SpatialParams<TypeTag, TTag::PorousMediumOnePNC> {
     using type = TwoPSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
 };
 #else
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+struct SpatialParams<TypeTag, TTag::PorousMediumOnePNC> {
     using type = OnePSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
 };
 #endif
 
 //! Use a model with constant tortuosity for the effective diffusivity
 template<class TypeTag>
-struct EffectiveDiffusivityModel<TypeTag, TTag::DarcyOnePNC>
+struct EffectiveDiffusivityModel<TypeTag, TTag::PorousMediumOnePNC>
 { using type = DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>; };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
 
 } //end namespace Dumux::Properties
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
index 716c9a8a926e2022591ab07cf68e9de5e0596e5f..2b9a725720b46d1f834ca791b2ce73852be2eeca 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -1,13 +1,16 @@
 dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_original
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=0)
 
 dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_a_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=1)
 
 dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_b_solution
                SOURCES main.cc
+               LABELS ffpm
                COMPILE_DEFINITIONS EXNUMBER=2)
 
 # add a symlink for each input file
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
index cfba40383d607ea9797221682cbe1fd7b2e1fe56..66dbbc3b192757f538d136af30c090d0496fc083 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
@@ -20,22 +20,25 @@
  * \file
  * \brief The free-flow sub problem
  */
-#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
-#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_SOL_HH
+#define DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_SOL_HH
 
-#if EXNUMBER >= 1
-#include <dumux/freeflow/rans/problem.hh>
-#include <dumux/freeflow/rans/boundarytypes.hh>
-#else
-#include <dumux/freeflow/navierstokes/staggered/problem.hh>
-#endif
 
 #include <dumux/common/timeloop.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
 #include <dumux/common/numeqvector.hh>
 #include <dumux/multidomain/boundary/stokesdarcy/couplingdata.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/freeflow/turbulencemodel.hh>
+#include <dumux/freeflow/turbulenceproperties.hh>
+#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
+#include <dumux/freeflow/rans/boundarytypes.hh>
+#else
+#include <dumux/freeflow/navierstokes/staggered/problem.hh>
 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
+#endif
 
 namespace Dumux {
 
@@ -54,6 +57,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
 #endif
 
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+    static constexpr auto dimWorld = GridView::dimensionworld;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -61,11 +65,11 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
 #if EXNUMBER >= 1
     using BoundaryTypes = Dumux::RANSBoundaryTypes<ModelTraits, ModelTraits::numEq()>;
 #else
-    using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag,
-    Properties::ModelTraits>::numEq()>;
+    using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
 #endif
     using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
@@ -85,8 +89,11 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
     static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
 
 public:
-    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+    FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                       std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "Freeflow"),
+    eps_(1e-6),
+    couplingManager_(couplingManager)
     {
         refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
         refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
@@ -95,17 +102,22 @@ public:
 
         diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
                                                                      getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
-    }
-
-
-   /*!
-     * \brief Return the sources within the domain.
-     *
-     * \param globalPos The global position
-     */
-    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
-    { return NumEqVector(0.0); }
+#if EXNUMBER >=1
+        FluidSystem::init();
+        Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
+        FluidState fluidState;
+        const auto phaseIdx = 0;
+        fluidState.setPressure(phaseIdx, refPressure_);
+        fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
+        fluidState.setMassFraction(phaseIdx, phaseIdx, refMoleFrac_);
+        Scalar density = FluidSystem::density(fluidState, phaseIdx);
+        Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, phaseIdx) / density;
+        Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+        turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(refVelocity_, diameter, kinematicViscosity);
+        dissipation_ = turbulenceProperties.dissipationRate(refVelocity_, diameter, kinematicViscosity);
+#endif
 
+    }
 
     /*!
      * \brief Specifies which kind of boundary condition should be
@@ -121,6 +133,20 @@ public:
 
         const auto& globalPos = scvf.center();
 
+#if EXNUMBER >=1
+        if(onRightBoundary_(globalPos))
+        {
+            values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
+            values.setOutflow(Indices::dissipationEqIdx);
+        }
+        else
+        {
+            // walls and inflow
+            values.setDirichlet(Indices::turbulentKineticEnergyIdx);
+            values.setDirichlet(Indices::dissipationIdx);
+        }
+#endif
+
         if (onLeftBoundary_(globalPos))
         {
             values.setDirichlet(Indices::velocityXIdx);
@@ -176,20 +202,70 @@ public:
             values.setCouplingNeumann(Indices::momentumYBalanceIdx);
             values.setBeaversJoseph(Indices::momentumXBalanceIdx);
         }
+
         return values;
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+     * \brief Evaluate the boundary conditions for a dirichlet values at the boundary.
      *
-     * \param element The element
-     * \param scvf The subcontrolvolume face
+     * \param element The finite element
+     * \param scvf the sub control volume face
+     * \note used for cell-centered discretization schemes
      */
-    PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+    PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
     {
-        PrimaryVariables values(0.0);
-        values = initialAtPos(pos);
+        const auto globalPos = scvf.ipGlobal();
+        PrimaryVariables values(initialAtPos(globalPos));
+        return values;
+    }
 
+    /*!
+     * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param scv The sub control volume
+     * \param pvIdx The primary variable index in the solution vector
+     */
+    template<class Element, class FVElementGeometry, class SubControlVolume>
+    bool isDirichletCell(const Element& element,
+                         const FVElementGeometry& fvGeometry,
+                         const SubControlVolume& scv,
+                         int pvIdx) const
+    {
+        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+        {
+            // For the komega model we set a fixed dissipation (omega) for all cells at the wall
+            for (const auto& scvf : scvfs(fvGeometry))
+                if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx)
+                    return true;
+        }
+        return false;
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for fixed values at cell centers
+     *
+     * \param element The finite element
+     * \param scv the sub control volume
+     * \note used for cell-centered discretization schemes
+     */
+    PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
+    {
+        const auto globalPos = scv.center();
+        PrimaryVariables values(initialAtPos(globalPos));
+        unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
+
+        if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+        {
+            // For the komega model we set a fixed value for the dissipation
+            const auto wallDistance = ParentType::wallDistance(elementIdx);
+            using std::pow;
+            values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
+                                                    / (ParentType::betaOmega() * wallDistance * wallDistance);
+            return values;
+        }
         return values;
     }
 
@@ -221,21 +297,15 @@ public:
         return values;
     }
 
-
     /*!
-     * \brief Set the coupling manager
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
      */
-    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
-    { couplingManager_ = cm; }
+    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+    { return NumEqVector(0.0); }
 
     /*!
-     * \brief Get the coupling manager
-     */
-    const CouplingManager& couplingManager() const
-    { return *couplingManager_; }
-
-
-   /*!
      * \brief Evaluate the initial value for a control volume.
      *
      * \param globalPos The global position
@@ -253,6 +323,25 @@ public:
         values[Indices::velocityXIdx] = refVelocity();
         values[Indices::temperatureIdx] = refTemperature();
 
+#if EXNUMBER >= 1
+        values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
+        values[Indices::dissipationIdx] = dissipation_;
+#endif
+
+#if EXNUMBER == 1
+        if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
+        {
+            values[Indices::turbulentKineticEnergyIdx] = 0.0;
+            values[Indices::dissipationIdx] = 0.0;
+        }
+#elif EXNUMBER >= 2
+        if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
+        {
+            values[Indices::turbulentKineticEnergyIdx] = 0.0;
+            values[Indices::dissipationIdx] = 0.0;
+        }
+#endif
+
 #if EXNUMBER >= 2
         if(onLowerBoundary_(globalPos))
             values[Indices::velocityXIdx] = 0.0;
@@ -280,7 +369,6 @@ public:
     const Scalar refTemperature() const
     { return refTemperature_; }
 
-
     void setTimeLoop(TimeLoopPtr timeLoop)
     { timeLoop_ = timeLoop; }
 
@@ -288,19 +376,25 @@ public:
      * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center()); }
 
     /*!
      * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
      */
     Scalar alphaBJ(const SubControlVolumeFace& scvf) const
-    {
-        return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
-    }
+    { return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
 
-    // \}
+    /*!
+     * \brief Set the coupling manager
+     */
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    /*!
+     * \brief Get the coupling manager
+     */
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
 
 private:
     bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -348,6 +442,8 @@ private:
     Scalar refPressure_;
     Scalar refMoleFrac_;
     Scalar refTemperature_;
+    Scalar turbulentKineticEnergy_;
+    Scalar dissipation_;
 
     TimeLoopPtr timeLoop_;
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc b/exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc
index c8974ad84a70a8b501b07f11dbbb73ae90c4daac..ad0470ad0a2755d83bee78e4ef1fc25a43fecac4 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc
@@ -62,51 +62,51 @@ int main(int argc, char** argv)
     Parameters::init(argc, argv);
 
     // Define the sub problem type tags
-    using StokesTypeTag = Properties::TTag::StokesZeroEq;
-    using DarcyTypeTag = Properties::TTag::DarcyTwoPTwoCNI;
+    using FreeflowTypeTag = Properties::TTag::FreeflowModel;
+    using PorousMediumTypeTag = Properties::TTag::PorousMediumModel;
 
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
-    using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
-    DarcyGridManager darcyGridManager;
-    darcyGridManager.init("Darcy"); // pass parameter group
+    using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+    PorousMediumGridManager porousMediumGridManager;
+    porousMediumGridManager.init("PorousMedium"); // pass parameter group
 
-    using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
-    StokesGridManager stokesGridManager;
-    stokesGridManager.init("Stokes"); // pass parameter group
+    using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+    FreeflowGridManager freeflowGridManager;
+    freeflowGridManager.init("Freeflow"); // pass parameter group
 
     // we compute on the leaf grid view
-    const auto& darcyGridView = darcyGridManager.grid().leafGridView();
-    const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+    const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+    const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
 
     // create the finite volume grid geometry
-    using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::GridGeometry>;
-    auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
-    using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>;
-    auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+    using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+    auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+    using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+    auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
 
-    using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+    using Traits = StaggeredMultiDomainTraits<FreeflowTypeTag, FreeflowTypeTag, PorousMediumTypeTag>;
 
     // the coupling manager
     using CouplingManager = StokesDarcyCouplingManager<Traits>;
-    auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+    auto couplingManager = std::make_shared<CouplingManager>(freeflowFvGridGeometry, porousMediumFvGridGeometry);
 
     // the indices
-    constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
-    constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
-    constexpr auto darcyIdx = CouplingManager::darcyIdx;
+    constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+    constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+    constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
 
     // the problem (initial and boundary conditions)
-    using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
-    auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
-    using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
-    auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+    using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+    auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+    using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+    auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
 
     // initialize the fluidsystem (tabulation)
-    GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+    GetPropType<FreeflowTypeTag, Properties::FluidSystem>::init();
 
     // get some time loop parameters
-    using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+    using Scalar = GetPropType<FreeflowTypeTag, Properties::Scalar>;
     const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
     const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -116,61 +116,61 @@ int main(int argc, char** argv)
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // set timeloop for the subproblems, needed for boundary value variations
-    stokesProblem->setTimeLoop(timeLoop);
-    darcyProblem->setTimeLoop(timeLoop);
+    freeflowProblem->setTimeLoop(timeLoop);
+    porousMediumProblem->setTimeLoop(timeLoop);
 
     // the solution vector
     Traits::SolutionVector sol;
-    sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
-    sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
-    sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+    sol[freeflowCellCenterIdx].resize(freeflowFvGridGeometry->numCellCenterDofs());
+    sol[freeflowFaceIdx].resize(freeflowFvGridGeometry->numFaceDofs());
+    sol[porousMediumIdx].resize(porousMediumFvGridGeometry->numDofs());
 
-    auto stokesSol = partial(sol, stokesFaceIdx, stokesCellCenterIdx);
+    auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
 
-    stokesProblem->applyInitialSolution(stokesSol);
-    darcyProblem->applyInitialSolution(sol[darcyIdx]);
+    freeflowProblem->applyInitialSolution(freeflowSol);
+    porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
 
     auto solOld = sol;
 
-    couplingManager->init(stokesProblem, darcyProblem, sol);
+    couplingManager->init(freeflowProblem, porousMediumProblem, sol);
     // TODO: update static wall properties
     // TODO: update dynamic wall properties
 #if EXNUMBER >= 1
-    stokesProblem->updateStaticWallProperties();
-    stokesProblem->updateDynamicWallProperties(stokesSol);
+    freeflowProblem->updateStaticWallProperties();
+    freeflowProblem->updateDynamicWallProperties(freeflowSol);
 #endif
 
     // the grid variables
-    using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
-    auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
-    stokesGridVariables->init(stokesSol);
-    using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
-    auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
-    darcyGridVariables->init(sol[darcyIdx]);
-
-    // initialize the vtk output module
-    const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
-    const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
-
-    StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
-    GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
-    stokesVtkWriter.write(0.0);
-
-    VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
-    using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
-    darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
-    GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
-    darcyVtkWriter.write(0.0);
+    using FreeflowGridVariables = GetPropType<FreeflowTypeTag, Properties::GridVariables>;
+    auto freeflowGridVariables = std::make_shared<FreeflowGridVariables>(freeflowProblem, freeflowFvGridGeometry);
+    freeflowGridVariables->init(freeflowSol);
+    using PorousMediumGridVariables = GetPropType<PorousMediumTypeTag, Properties::GridVariables>;
+    auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
+    porousMediumGridVariables->init(sol[porousMediumIdx]);
+
+    // intialize the vtk output module
+    const auto freeflowName = getParam<std::string>("Problem.Name") + "_" + freeflowProblem->name();
+    const auto porousMediumName = getParam<std::string>("Problem.Name") + "_" + porousMediumProblem->name();
+
+    StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
+    GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
+    freeflowVtkWriter.write(0.0);
+
+    VtkOutputModule<PorousMediumGridVariables, GetPropType<PorousMediumTypeTag, Properties::SolutionVector>> porousMediumVtkWriter(*porousMediumGridVariables, sol[porousMediumIdx], porousMediumName);
+    using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
+    porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
+    GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
+    porousMediumVtkWriter.write(0.0);
 
     // the assembler with time loop for instationary problem
     using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
-    auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
-                                                 std::make_tuple(stokesFvGridGeometry->faceFVGridGeometryPtr(),
-                                                                 stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
-                                                                 darcyFvGridGeometry),
-                                                 std::make_tuple(stokesGridVariables->faceGridVariablesPtr(),
-                                                                 stokesGridVariables->cellCenterGridVariablesPtr(),
-                                                                 darcyGridVariables),
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(freeflowProblem, freeflowProblem, porousMediumProblem),
+                                                 std::make_tuple(freeflowFvGridGeometry->faceFVGridGeometryPtr(),
+                                                                 freeflowFvGridGeometry->cellCenterFVGridGeometryPtr(),
+                                                                 porousMediumFvGridGeometry),
+                                                 std::make_tuple(freeflowGridVariables->faceGridVariablesPtr(),
+                                                                 freeflowGridVariables->cellCenterGridVariablesPtr(),
+                                                                 porousMediumGridVariables),
                                                  couplingManager,
                                                  timeLoop, solOld);
 
@@ -193,22 +193,22 @@ int main(int argc, char** argv)
 
 #if EXNUMBER >= 1
         // TODO: update dynamic wall properties
-        stokesProblem->updateDynamicWallProperties(stokesSol);
+        freeflowProblem->updateDynamicWallProperties(freeflowSol);
 #endif
 
-        // post time step treatment of Darcy problem
-        darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+        // post time step treatment of PorousMedium problem
+        porousMediumProblem->postTimeStep(sol[porousMediumIdx], *porousMediumGridVariables, timeLoop->timeStepSize());
 
         // advance grid variables to the next time step
-        stokesGridVariables->advanceTimeStep();
-        darcyGridVariables->advanceTimeStep();
+        freeflowGridVariables->advanceTimeStep();
+        porousMediumGridVariables->advanceTimeStep();
 
         // advance to the time loop to the next step
         timeLoop->advanceTimeStep();
 
         // write vtk output
-        stokesVtkWriter.write(timeLoop->time());
-        darcyVtkWriter.write(timeLoop->time());
+        freeflowVtkWriter.write(timeLoop->time());
+        porousMediumVtkWriter.write(timeLoop->time());
 
         // report statistics of this time step
         timeLoop->reportTimeStep();
@@ -218,8 +218,8 @@ int main(int argc, char** argv)
 
     } while (!timeLoop->finished());
 
-    timeLoop->finalize(stokesGridView.comm());
-    timeLoop->finalize(darcyGridView.comm());
+    timeLoop->finalize(freeflowGridView.comm());
+    timeLoop->finalize(porousMediumGridView.comm());
 
     Parameters::print();
 
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
index 5ab29c8c810135af206cd4bf210e223cb28eaab7..0769cc10d4f3977f7a3b0bd6c25bcac56c122cc9 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
@@ -3,7 +3,7 @@ DtInitial =  1e-1 # [s]
 MaxTimeStepSize = 43200 # [s] (12 hours)
 TEnd = 864000 # [s] (6 days)
 
-[Stokes.Grid]
+[Freeflow.Grid]
 Positions0 = 0.0 0.25
 Positions1 = 0.25 0.5
 Grading0 = 1.0
@@ -12,7 +12,7 @@ Cells0 = 15
 Cells1 = 20
 Verbosity = true
 
-[Darcy.Grid]
+[PorousMedium.Grid]
 Positions0 = 0.0 0.25
 Positions1 = 0.0 0.25
 Cells0 = 15
@@ -21,16 +21,16 @@ Grading0 = 1.0
 Grading1 = 1.0
 Verbosity = true
 
-[Stokes.Problem]
-Name = stokes
+[Freeflow.Problem]
+Name = freeflow
 RefVelocity = 3.5 # [m/s]
 RefPressure = 1e5 # [Pa]
 refMoleFrac = 0 # [-]
 RefTemperature = 298.15 # [K]
 EnableInertiaTerms = true
 
-[Darcy.Problem]
-Name = darcy
+[PorousMedium.Problem]
+Name = porousmedium
 Pressure = 1e5
 Saturation = 0.5 # initial Sw
 Temperature = 298.15 # [K]
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
index 456f9b14fac589d5b233dfb87c073f5ad3a054d5..f7fd1fd835a6f3cebae701e5f09dc8ff578618c2 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/porousmediumsubproblem.hh
@@ -21,8 +21,8 @@
  *
  * \brief The porous medium sub problem
  */
-#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
-#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#ifndef DUMUX_POROUSMEDIUMFLOW2P2C_SUBPROBLEM_HH
+#define DUMUX_POROUSMEDIUMFLOW2P2C_SUBPROBLEM_HH
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/boundarytypes.hh>
@@ -77,8 +77,10 @@ class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
 
 public:
     PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
-                   std::shared_ptr<CouplingManager> couplingManager)
-    : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+                           std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry, "PorousMedium"),
+    eps_(1e-7),
+    couplingManager_(couplingManager)
     {
         pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
         initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
@@ -123,12 +125,12 @@ public:
     }
 
     /*!
-      * \brief Specifies which kind of boundary condition should be
-      *        used for which equation on a given boundary control volume.
-      *
-      * \param element The element
-      * \param scvf The boundary sub control volume face
-      */
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param element The element
+     * \param scvf The boundary sub control volume face
+     */
     BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
     {
         BoundaryTypes values;
@@ -207,8 +209,6 @@ public:
                        const SubControlVolume &scv) const
     { return NumEqVector(0.0); }
 
-    // \}
-
     /*!
      * \brief Evaluate the initial value for a control volume.
      *
@@ -227,8 +227,6 @@ public:
         return values;
     }
 
-    // \}
-
     /*!
      * \brief Set the coupling manager
      */
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
index d0108365644dbf0a584544168a6c8976fde7e8a4..3b6b962a2ced962347b31f66e892739f3bafa869 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
@@ -43,7 +43,7 @@
 #include <dumux/discretization/staggered/freeflow/properties.hh>
 
 #if EXNUMBER >= 1
-#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/compositional/sstncmodel.hh>
 #else
 #include <dumux/freeflow/compositional/navierstokesncmodel.hh>
 #endif
@@ -54,78 +54,78 @@ namespace Dumux::Properties {
 
 // Create new type tags
 namespace TTag {
-struct DarcyTwoPTwoCNI { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+struct PorousMediumModel { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
 #if EXNUMBER >= 1
-struct StokesZeroEq { using InheritsFrom = std::tuple<ZeroEqNCNI, StaggeredFreeFlowModel>; };
+struct FreeflowModel { using InheritsFrom = std::tuple<SSTNCNI, StaggeredFreeFlowModel>; };
 #else
-struct StokesZeroEq { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+struct FreeflowModel { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
 #endif
 } // end namespace TTag
 
 // Set the coupling manager
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::StokesZeroEq>
+struct CouplingManager<TypeTag, TTag::FreeflowModel>
 {
-    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCNI>;
+    using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumModel>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct CouplingManager<TypeTag, TTag::PorousMediumModel>
 {
-    using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesZeroEq, Properties::TTag::StokesZeroEq, TypeTag>;
+    using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowModel, Properties::TTag::FreeflowModel, TypeTag>;
     using type = Dumux::StokesDarcyCouplingManager<Traits>;
 };
 
 // Set the problem property
 template<class TypeTag>
-struct Problem<TypeTag, TTag::StokesZeroEq> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+struct Problem<TypeTag, TTag::FreeflowModel> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
 template<class TypeTag>
-struct Problem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
+struct Problem<TypeTag, TTag::PorousMediumModel> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
 
 // Set the grid type
 template<class TypeTag>
-struct Grid<TypeTag, TTag::StokesZeroEq> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::FreeflowModel> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 template<class TypeTag>
-struct Grid<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+struct Grid<TypeTag, TTag::PorousMediumModel> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
 
 // The fluid system
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::StokesZeroEq>
+struct FluidSystem<TypeTag, TTag::FreeflowModel>
 {
   using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
   static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
   using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
 };
 template<class TypeTag>
-struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
+struct FluidSystem<TypeTag, TTag::PorousMediumModel> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
 
 // The gas component balance (air) is replaced by the total mass balance
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::PorousMediumModel> { static constexpr int value = 3; };
 template<class TypeTag>
-struct ReplaceCompEqIdx<TypeTag, TTag::StokesZeroEq> { static constexpr int value = 3; };
+struct ReplaceCompEqIdx<TypeTag, TTag::FreeflowModel> { static constexpr int value = 3; };
 
 // Use formulation based on mass fractions
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr bool value = true; };
+struct UseMoles<TypeTag, TTag::PorousMediumModel> { static constexpr bool value = true; };
 
 //! Set the default formulation to pw-Sn: This can be over written in the problem.
 template<class TypeTag>
-struct Formulation<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct Formulation<TypeTag, TTag::PorousMediumModel>
 { static constexpr auto value = TwoPFormulation::p1s0; };
 
 template<class TypeTag>
-struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCNI>
+struct SpatialParams<TypeTag, TTag::PorousMediumModel>
 { using type = TwoPSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>; };
 
 template<class TypeTag>
-struct EnableGridGeometryCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridGeometryCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 template<class TypeTag>
-struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowModel> { static constexpr bool value = true; };
 
 } //end namespace Dumux::Properties