diff --git a/doc/doxygen/modules.txt b/doc/doxygen/modules.txt
index 2414dbe80966720e2dfc061fd074d3d0049d9651..52cf17f1410ac952504352c24cb0a64bba33cb99 100644
--- a/doc/doxygen/modules.txt
+++ b/doc/doxygen/modules.txt
@@ -207,6 +207,42 @@
          *
          * \copydetails Dumux::StokesncniModel
          */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqModel ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqModel
+         */
+         /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqncModel N-component ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqncModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqncniModel Non-isothermal N-component ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqncniModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqModel ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqModel
+         */
+         /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqncModel N-component ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqncModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup BoxZeroEqncniModel Non-isothermal N-component ZeroEq
+         *
+         * \copydetails Dumux::ZeroEqncniModel
+         */
         /*!
          * \ingroup ImplicitModels
          * \defgroup ElasticBoxModel Linear elastic 
@@ -233,16 +269,42 @@
          * \ingroup ImplicitModels
          * \defgroup TwoPTwoCStokesTwoCModel Two-component, Stokes-Darcy
          *
+         * \copydetails Dumux::TwoCStokesTwoPTwoCLocalOperator
+         * <br><br><br>
          * \copydetails Dumux::TwoPTwoCModel
+         * <br>
          * \copydetails Dumux::StokesncModel
          */
         /*!
          * \ingroup ImplicitModels
          * \defgroup TwoPTwoCNIStokesTwoCNIModel Non-isothermal, two-component, Stokes-Darcy
          *
+         * \copydetails Dumux::TwoCNIStokesTwoPTwoCNILocalOperator
+         * <br><br><br>
          * \copydetails Dumux::TwoPTwoCNIModel
+         * <br>
          * \copydetails Dumux::StokesncniModel
          */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup TwoPTwoCZeroEqTwoCModel Two-component, ZeroEq turbulence-Darcy
+         *
+         * \copydetails Dumux::TwoCStokesTwoPTwoCLocalOperator
+         * <br><br><br>
+         * \copydetails Dumux::TwoPTwoCModel
+         * <br>
+         * \copydetails Dumux::ZeroEqncModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup TwoPTwoCNIZeroEqTwoCNIModel Non-isothermal, two-component, ZeroEq turbulence-Darcy
+         *
+         * \copydetails Dumux::TwoCNIStokesTwoPTwoCNILocalOperator
+         * <br><br><br>
+         * \copydetails Dumux::TwoPTwoCNIModel
+         * <br>
+         * \copydetails Dumux::ZeroEqncniModel
+         */
     /*!
      * \ingroup ImplicitModel
      * \defgroup ImplicitBaseProblems Base Problems
diff --git a/doc/handbook/CMakeLists.txt b/doc/handbook/CMakeLists.txt
index 4cb77b9e62d63ab2f8c836b205744b66de68c91c..bbd7bcb73175a0dde9a319c4bed2fb2978c58e18 100644
--- a/doc/handbook/CMakeLists.txt
+++ b/doc/handbook/CMakeLists.txt
@@ -20,6 +20,10 @@ set(TEX_INPUTS
   tutorial-coupled.tex
   tutorial-decoupled.tex
   tutorial.tex
+  ModelDescriptions/2cstokes2p2cmodel.tex
+  ModelDescriptions/2cnistokes2p2cnimodel.tex
+  ModelDescriptions/2czeroeq2p2cmodel.tex
+  ModelDescriptions/2cnizeroeq2p2cnimodel.tex
   ModelDescriptions/1p2cimplicitmodel.tex
   ModelDescriptions/1pdecoupledmodel.tex
   ModelDescriptions/1pimplicitmodel.tex
@@ -42,6 +46,9 @@ set(TEX_INPUTS
   ModelDescriptions/stokesimplicitmodel.tex
   ModelDescriptions/stokesncimplicitmodel.tex
   ModelDescriptions/stokesncniimplicitmodel.tex
+  ModelDescriptions/zeroeqimplicitmodel.tex
+  ModelDescriptions/zeroeqncimplicitmodel.tex
+  ModelDescriptions/zeroeqncniimplicitmodel.tex
   ../../tutorial/tutorial_coupled.cc
   ../../tutorial/tutorial_coupled.input
   ../../tutorial/tutorialproblem_coupled.hh
diff --git a/doc/handbook/ModelDescriptions/2cnistokes2p2cnimodel.tex b/doc/handbook/ModelDescriptions/2cnistokes2p2cnimodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..6505b25ab059d6586d02a158cdffda8d45a09063
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/2cnistokes2p2cnimodel.tex
@@ -0,0 +1,24 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements the coupling between a free-\/flow model and a porous-\/medium flow model under non-\/isothermal conditions. Here the coupling conditions for the individual balance are presented\-:
+
+The total mass balance equation\-: \[ \left[ \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = -\left[ \left( \varrho_\textrm{g} \boldsymbol{v}_\textrm{g} + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} \right) \cdot \boldsymbol{n} \right]^\textrm{pm} \] in which $n$ represents a vector normal to the interface pointing outside of the specified subdomain.
+
+The momentum balance (tangential), which corresponds to the Beavers-\/\-Jospeh Saffman condition\-: \[ \left[ \left( {\boldsymbol{v}}_\textrm{g} + \frac{\sqrt{\left(\boldsymbol{K} \boldsymbol{t}_i \right) \cdot \boldsymbol{t}_i}} {\alpha_\textrm{BJ} \mu_\textrm{g}} \boldsymbol{{\tau}}_\textrm{t} \boldsymbol{n} \right) \cdot \boldsymbol{t}_i \right]^\textrm{ff} = 0 \] with $ \boldsymbol{{\tau}_\textrm{t}} = \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right] \nabla \left( \boldsymbol{v}_\textrm{g} + \boldsymbol{v}_\textrm{g}^\intercal \right) $ in which the eddy viscosity $ \mu_\textrm{g,t} = 0 $ for the Stokes equation.
+
+The momentum balance (normal)\-: \[ \left[ \left( \left\lbrace \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} {\boldsymbol{v}}_\textrm{g}^\intercal - \boldsymbol{{\tau}}_\textrm{t} + {p}_\textrm{g} \boldsymbol{I} \right\rbrace \boldsymbol{n} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = p_\textrm{g}^\textrm{pm} \]
+
+The component mass balance equation (continuity of fluxes)\-: \[ \left[ \left( \varrho_\textrm{g} {X}^\kappa_\textrm{g} {\boldsymbol{v}}_\textrm{g} - {\boldsymbol{j}}^\kappa_\textrm{g,ff,t,diff} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = -\left[ \left( \varrho_\textrm{g} X^\kappa_\textrm{g} \boldsymbol{v}_\textrm{g} - \boldsymbol{j}^\kappa_\textrm{g,pm,diff} + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} X^\kappa_\textrm{l} - \boldsymbol{j}^\kappa_\textrm{l,pm,diff} \right) \cdot \boldsymbol{n} \right]^\textrm{pm} = 0 \] in which the diffusive fluxes $ j_\textrm{diff} $ are the diffusive fluxes as they are implemented in the individual subdomain models.
+
+The component mass balance equation (continuity of mass/ mole fractions)\-: \[ \left[ {X}^{\kappa}_\textrm{g} \right]^\textrm{ff} = \left[ X^{\kappa}_\textrm{g} \right]^\textrm{pm} \]
+
+The energy balance equation (continuity of fluxes)\-: \[ \left[ \left( \varrho_\textrm{g} {h}_\textrm{g} {\boldsymbol{v}}_\textrm{g} - {h}^\textrm{a}_\textrm{g} {\boldsymbol{j}}^\textrm{a}_\textrm{g,ff,t,diff} - {h}^\textrm{w}_\textrm{g} {\boldsymbol{j}}^\textrm{w}_\textrm{g,ff,t,diff} - \left( \lambda_\textrm{g} + \lambda_\textrm{g,t} \right) \nabla {T} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = -\left[ \left( \varrho_\textrm{g} h_\textrm{g} \boldsymbol{v}_\textrm{g} + \varrho_\textrm{l} h_\textrm{l} \boldsymbol{v}_\textrm{l} - \lambda_\textrm{pm} \nabla T \right) \cdot \boldsymbol{n} \right]^\textrm{pm} \]
+
+The energy balance equation (continuity of temperature)\-: \[ \left[ {T} \right]^\textrm{ff} = \left[ T \right]^\textrm{pm} \]
+
+This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
+
diff --git a/doc/handbook/ModelDescriptions/2cnizeroeq2p2cnimodel.tex b/doc/handbook/ModelDescriptions/2cnizeroeq2p2cnimodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..81d55e0a120bd1d88ef3550e96f7a8707e185ac7
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/2cnizeroeq2p2cnimodel.tex
@@ -0,0 +1,6 @@
+
+This model implements the same coupling conditions as the model
+concepts, which couples Darcy and Stokes flow presented \hyperref[sc_2cnistokes2p2cni]{above}.
+The only difference is that the eddy coefficients $\mu_\textrm{g,t}$,
+$D_\textrm{g,t}$, and $\lambda_\textrm{g,t}$ are not necessarily non-zero
+at the coupling interface.
diff --git a/doc/handbook/ModelDescriptions/2cstokes2p2cmodel.tex b/doc/handbook/ModelDescriptions/2cstokes2p2cmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..55e263dbb97b06dc8a5c06066a78e902ce65ec6e
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/2cstokes2p2cmodel.tex
@@ -0,0 +1,20 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements the coupling between a free-\/flow model and a porous-\/medium flow model under isothermal conditions. Here the coupling conditions for the individual balance are presented\-:
+
+The total mass balance equation\-: \[ \left[ \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = -\left[ \left( \varrho_\textrm{g} \boldsymbol{v}_\textrm{g} + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} \right) \cdot \boldsymbol{n} \right]^\textrm{pm} \] in which $n$ represents a vector normal to the interface pointing outside of the specified subdomain.
+
+The momentum balance (tangential), which corresponds to the Beavers-\/\-Jospeh Saffman condition\-: \[ \left[ \left( {\boldsymbol{v}}_\textrm{g} + \frac{\sqrt{\left(\boldsymbol{K} \boldsymbol{t}_i \right) \cdot \boldsymbol{t}_i}} {\alpha_\textrm{BJ} \mu_\textrm{g}} \boldsymbol{{\tau}}_\textrm{t} \boldsymbol{n} \right) \cdot \boldsymbol{t}_i \right]^\textrm{ff} = 0 \] with $ \boldsymbol{{\tau}_\textrm{t}} = \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right] \nabla \left( \boldsymbol{v}_\textrm{g} + \boldsymbol{v}_\textrm{g}^\intercal \right) $ in which the eddy viscosity $ \mu_\textrm{g,t} = 0 $ for the Stokes equation.
+
+The momentum balance (normal)\-: \[ \left[ \left( \left\lbrace \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} {\boldsymbol{v}}_\textrm{g}^\intercal - \boldsymbol{{\tau}}_\textrm{t} + {p}_\textrm{g} \boldsymbol{I} \right\rbrace \boldsymbol{n} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = p_\textrm{g}^\textrm{pm} \]
+
+The component mass balance equation (continuity of fluxes)\-: \[ \left[ \left( \varrho_\textrm{g} {X}^\kappa_\textrm{g} {\boldsymbol{v}}_\textrm{g} - {\boldsymbol{j}}^\kappa_\textrm{g,ff,t,diff} \right) \cdot \boldsymbol{n} \right]^\textrm{ff} = -\left[ \left( \varrho_\textrm{g} X^\kappa_\textrm{g} \boldsymbol{v}_\textrm{g} - \boldsymbol{j}^\kappa_\textrm{g,pm,diff} + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} X^\kappa_\textrm{l} - \boldsymbol{j}^\kappa_\textrm{l,pm,diff} \right) \cdot \boldsymbol{n} \right]^\textrm{pm} = 0 \] in which the diffusive fluxes $ j_\textrm{diff} $ are the diffusive fluxes as they are implemented in the individual subdomain models.
+
+The component mass balance equation (continuity of mass/ mole fractions)\-: \[ \left[ {X}^{\kappa}_\textrm{g} \right]^\textrm{ff} = \left[ X^{\kappa}_\textrm{g} \right]^\textrm{pm} \]
+
+This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
+
diff --git a/doc/handbook/ModelDescriptions/2czeroeq2p2cmodel.tex b/doc/handbook/ModelDescriptions/2czeroeq2p2cmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..172f325fd863fa6a6bb064bda747939184cccbd4
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/2czeroeq2p2cmodel.tex
@@ -0,0 +1,6 @@
+
+This model implements the same coupling conditions as the model
+concepts, which couples Darcy and Stokes flow presented \hyperref[sc_2cstokes2p2c]{above}.
+The only difference is that the eddy coefficients $\mu_\textrm{g,t}$
+and $D_\textrm{g,t}$ are not necessarily non-zero at the coupling
+interface.
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
index 5bf2a4afac7d31bbe6163aeafaee7b7d2de3f7a1..42f7be883ccbf37590ee9690d06ed61ef43e31df 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
@@ -8,7 +8,7 @@ This model solves equations of the form \[ \phi \left( \rho_w \frac{\partial S_w
 
 \[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_w + f_n \textbf{grad}\, p_c + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q, \]
 
-a non-\/wetting ( $ n $) phase pressure yields \[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_n - f_w \textbf{grad}\, p_c + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q, \] and a global pressure leads to \[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_{global} + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q. \] Here, $ p_\alpha $ is a phase pressure, $ p_ {global} $ the global pressure of a classical fractional flow formulation (see e.\-g. P. Binning and M. A. Celia, ''Practical implementation of the fractional flow approach to multi-\/phase flow simulation'', Advances in water resources, vol. 22, no. 5, pp. 461-\/478, 1999.), $ p_c = p_n - p_w $ is the capillary pressure, $ \boldsymbol K $ the absolute permeability, $ \lambda = \lambda_w + \lambda_n $ the total mobility depending on the saturation ( $ \lambda_\alpha = k_{r_\alpha} / \mu_\alpha $), $ f_\alpha = \lambda_\alpha / \lambda $ the fractional flow function of a phase, $ \rho_\alpha $ a phase density, $ g $ the gravity constant and $ q $ the source term.
+a non-\/wetting ( $ n $) phase pressure yields \[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_n - f_w \textbf{grad}\, p_c + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q, \] and a global pressure leads to \[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_{global} + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q. \] Here, $ p_\alpha $ is a phase pressure, $ p_ {global} $ the global pressure of a classical fractional flow formulation (see e.\-g. P. Binning and M. A. Celia, \textquotesingle{}\textquotesingle{}Practical implementation of the fractional flow approach to multi-\/phase flow simulation\textquotesingle{}\textquotesingle{}, Advances in water resources, vol. 22, no. 5, pp. 461-\/478, 1999.), $ p_c = p_n - p_w $ is the capillary pressure, $ \boldsymbol K $ the absolute permeability, $ \lambda = \lambda_w + \lambda_n $ the total mobility depending on the saturation ( $ \lambda_\alpha = k_{r_\alpha} / \mu_\alpha $), $ f_\alpha = \lambda_\alpha / \lambda $ the fractional flow function of a phase, $ \rho_\alpha $ a phase density, $ g $ the gravity constant and $ q $ the source term.
 
 For all cases, $ p = p_D $ on $ \Gamma_{Dirichlet} $, and $ \boldsymbol v_{total} \cdot \boldsymbol n = q_N $ on $ \Gamma_{Neumann} $.
 
@@ -18,7 +18,7 @@ In the I\-M\-P\-E\-S models the default setting is\-:
 
 
 \begin{itemize}
-\item formulation\-: $ p_w-S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00099_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
+\item formulation\-: $ p_w-S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00095_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
 \item compressibility\-: disabled (Property\-: {\itshape Enable\-Compressibility} set to {\itshape false})
 \end{itemize}
 
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
index b50f0704456ba30511d32ef7c2071512a2cc479c..8bf9b712ef2ceafe86eaaadc9ffc8a18c031a691 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
@@ -16,7 +16,7 @@ The total velocity formulation is only implemented for incompressible fluids and
 
 In the I\-M\-P\-E\-S models the default setting is\-:
 
-formulation\-: $ p_w $ -\/ $ S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00099_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
+formulation\-: $ p_w $ -\/ $ S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00095_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
 
 compressibility\-: disabled (Property\-: {\itshape Enable\-Compressibility} set to {\itshape false})
 
diff --git a/doc/handbook/ModelDescriptions/co2implicitmodel.tex b/doc/handbook/ModelDescriptions/co2implicitmodel.tex
index fb91c3715813a1d2a6ff10e045877d835bec6815..7812c03ac83cef8640a943de0f9c2a14a2fc8ca1 100644
--- a/doc/handbook/ModelDescriptions/co2implicitmodel.tex
+++ b/doc/handbook/ModelDescriptions/co2implicitmodel.tex
@@ -4,5 +4,5 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-See \hyperlink{a00615}{Two\-P\-Two\-C\-Model} for reference to the equations used. The \hyperlink{a00075}{C\-O2} model is derived from the 2p2c model. In the \hyperlink{a00075}{C\-O2} model the phase switch criterion is different from the 2p2c model. The phase switch occurs when the equilibrium concentration of a component in a phase is exceeded, instead of the sum of the components in the virtual phase (the phase which is not present) being greater that unity as done in the 2p2c model. The \hyperlink{a00082}{C\-O2\-Volume\-Variables} do not use a constraint solver for calculating the mole fractions as is the case in the 2p2c model. Instead mole fractions are calculated in the Fluid\-System with a given temperature, pressurem and salinity. The model is able to use either mole or mass fractions. The property use\-Moles can be set to either true or false in the problem file. Make sure that the according units are used in the problem setup. use\-Moles is set to false by default.
+See \hyperlink{a00633}{Two\-P\-Two\-C\-Model} for reference to the equations used. The \hyperlink{a00074}{C\-O2} model is derived from the 2p2c model. In the \hyperlink{a00074}{C\-O2} model the phase switch criterion is different from the 2p2c model. The phase switch occurs when the equilibrium concentration of a component in a phase is exceeded, instead of the sum of the components in the virtual phase (the phase which is not present) being greater that unity as done in the 2p2c model. The \hyperlink{a00078}{C\-O2\-Volume\-Variables} do not use a constraint solver for calculating the mole fractions as is the case in the 2p2c model. Instead mole fractions are calculated in the Fluid\-System with a given temperature, pressurem and salinity. The model is able to use either mole or mass fractions. The property use\-Moles can be set to either true or false in the problem file. Make sure that the according units are used in the problem setup. use\-Moles is set to false by default.
 
diff --git a/doc/handbook/ModelDescriptions/el1p2cmodel.tex b/doc/handbook/ModelDescriptions/el1p2cmodel.tex
index 23bdbac47efe11a0cac2b609b643107e9dbbe119..156dd70f5792e44509b4459fb418bf0bf0d91ba2 100644
--- a/doc/handbook/ModelDescriptions/el1p2cmodel.tex
+++ b/doc/handbook/ModelDescriptions/el1p2cmodel.tex
@@ -4,9 +4,9 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements a one-\/phase flow of an incompressible fluid, that consists of two components. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
+This model implements a one-\/phase flow of an incompressible fluid, that consists of two components. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke\textquotesingle{}s law.
 
-As an equation for the conservation of momentum within the fluid phase Darcy's approach is used\-: \[ v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right) \]
+As an equation for the conservation of momentum within the fluid phase Darcy\textquotesingle{}s approach is used\-: \[ v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right) \]
 
 Gravity can be enabled or disabled via the property system. By inserting this into the volume balance of the solid-\/fluid mixture, one gets \[ \frac{\partial \text{div} \textbf{u}}{\partial t} - \text{div} \left\{ \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right)\right\} = q \;, \]
 
diff --git a/doc/handbook/ModelDescriptions/el2pmodel.tex b/doc/handbook/ModelDescriptions/el2pmodel.tex
index d95cdf204c4c11b30b2e60945220dc9e624cc939..f60dd9a337e0619539d19bee9247c5746d4c6519 100644
--- a/doc/handbook/ModelDescriptions/el2pmodel.tex
+++ b/doc/handbook/ModelDescriptions/el2pmodel.tex
@@ -4,9 +4,9 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements a two-\/phase flow of compressible immiscible fluids $\alpha \in \{ w, n \}$. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
+This model implements a two-\/phase flow of compressible immiscible fluids $\alpha \in \{ w, n \}$. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke\textquotesingle{}s law.
 
-As an equation for the conservation of momentum within the fluid phases the standard multiphase Darcy's approach is used\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right) \]
+As an equation for the conservation of momentum within the fluid phases the standard multiphase Darcy\textquotesingle{}s approach is used\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right) \]
 
 Gravity can be enabled or disabled via the property system. By inserting this into the continuity equation, one gets \[ \frac{\partial \phi_{eff} \varrho_\alpha S_\alpha}{\partial t} - \text{div} \left\{ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}_\text{eff} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right) - \phi_{eff} \varrho_\alpha S_\alpha \frac{\partial \mathbf{u}}{\partial t} \right\} - q_\alpha = 0 \;, \]
 
diff --git a/doc/handbook/ModelDescriptions/elasticmodel.tex b/doc/handbook/ModelDescriptions/elasticmodel.tex
index 79fa14958672f889b8afcb676cefbd59caf2fab5..5c22a3c674581a91581a7915d4d4f1cc1125f5d4 100644
--- a/doc/handbook/ModelDescriptions/elasticmodel.tex
+++ b/doc/handbook/ModelDescriptions/elasticmodel.tex
@@ -4,7 +4,7 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements a linear elastic solid using Hooke's law as stress-\/strain relation and a quasi-\/stationary momentum balance equation\-: \[ \boldsymbol{\sigma} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \boldsymbol{I}. \]
+This model implements a linear elastic solid using Hooke\textquotesingle{}s law as stress-\/strain relation and a quasi-\/stationary momentum balance equation\-: \[ \boldsymbol{\sigma} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \boldsymbol{I}. \]
 
 with the strain tensor $\boldsymbol{\epsilon}$ as a function of the solid displacement gradient $\textbf{grad} \boldsymbol{u}$\-: \[ \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \boldsymbol{u} + \textbf{grad}^T \boldsymbol{u}). \]
 
diff --git a/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex b/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
index 5797e89418a53ce327518bbc5e72c99dba61a494..205aa3d9a41b9c3ca9a590609f0dc59089852330 100644
--- a/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
+++ b/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
@@ -6,7 +6,7 @@
 
 This model implements a $M$-\/phase flow of a fluid mixture composed of $N$ chemical species. The phases are denoted by lower index $\alpha \in \{ 1, \dots, M \}$. All fluid phases are mixtures of $N \geq M - 1$ chemical species which are denoted by the upper index $\kappa \in \{ 1, \dots, N \} $.
 
-The momentum approximation can be selected via \char`\"{}\-Base\-Flux\-Variables\char`\"{}\-: Darcy (\hyperlink{a00278}{Implicit\-Darcy\-Flux\-Variables}) and Forchheimer (\hyperlink{a00279}{Implicit\-Forchheimer\-Flux\-Variables}) relations are available for all Box models.
+The momentum approximation can be selected via \char`\"{}\-Base\-Flux\-Variables\char`\"{}\-: Darcy (\hyperlink{a00280}{Implicit\-Darcy\-Flux\-Variables}) and Forchheimer (\hyperlink{a00281}{Implicit\-Forchheimer\-Flux\-Variables}) relations are available for all Box models.
 
 By inserting this into the equations for the conservation of the mass of each component, one gets one mass-\/continuity equation for each component $\kappa$ \[ \sum_{\kappa} \left( \phi \frac{\partial \left(\varrho_\alpha x_\alpha^\kappa S_\alpha\right)}{\partial t} + \mathrm{div}\; \left\{ v_\alpha \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa \right\} \right) = q^\kappa \] with $\overline M_\alpha$ being the average molar mass of the phase $\alpha$\-: \[ \overline M_\alpha = \sum_\kappa M^\kappa \; x_\alpha^\kappa \]
 
@@ -20,7 +20,7 @@ These three equations constitute a non-\/linear complementarity problem, which c
 
 Several non-\/linear complementarity functions have been suggested, e.\-g. the Fischer-\/\-Burmeister function \[ \Phi(a,b) = a + b - \sqrt{a^2 + b^2} \;. \] This model uses \[ \Phi(a,b) = \min \{a, b \}\;, \] because of its piecewise linearity.
 
-These equations are then discretized using a fully-\/implicit vertex centered finite volume scheme (often known as 'box'-\/scheme) for spatial discretization and the implicit Euler method as temporal discretization.
+These equations are then discretized using a fully-\/implicit vertex centered finite volume scheme (often known as \textquotesingle{}box\textquotesingle{}-\/scheme) for spatial discretization and the implicit Euler method as temporal discretization.
 
 The model assumes local thermodynamic equilibrium and uses the following primary variables\-:
 \begin{itemize}
diff --git a/doc/handbook/ModelDescriptions/richardsimplicitmodel.tex b/doc/handbook/ModelDescriptions/richardsimplicitmodel.tex
index 526c0cda2161ea42d0915787868056efff1914a3..1d4cc399a02ef452a566e3aca69adaddab352faa 100644
--- a/doc/handbook/ModelDescriptions/richardsimplicitmodel.tex
+++ b/doc/handbook/ModelDescriptions/richardsimplicitmodel.tex
@@ -4,9 +4,9 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-In the unsaturated zone, Richards' equation \[ \frac{\partial\;\phi S_w \varrho_w}{\partial t} - \text{div} \left\lbrace \varrho_w \frac{k_{rw}}{\mu_w} \; \mathbf{K} \; \left( \text{\textbf{grad}} p_w - \varrho_w \textbf{g} \right) \right\rbrace = q_w, \] is frequently used to approximate the water distribution above the groundwater level.
+In the unsaturated zone, Richards\textquotesingle{} equation \[ \frac{\partial\;\phi S_w \varrho_w}{\partial t} - \text{div} \left\lbrace \varrho_w \frac{k_{rw}}{\mu_w} \; \mathbf{K} \; \left( \text{\textbf{grad}} p_w - \varrho_w \textbf{g} \right) \right\rbrace = q_w, \] is frequently used to approximate the water distribution above the groundwater level.
 
-It can be derived from the two-\/phase equations, i.\-e. \[ \phi\frac{\partial S_\alpha \varrho_\alpha}{\partial t} - \text{div} \left\lbrace \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}\; \mathbf{K} \; \left( \text{\textbf{grad}} p_\alpha - \varrho_\alpha \textbf{g} \right) \right\rbrace = q_\alpha, \] where $\alpha \in \{w, n\}$ is the fluid phase, $\kappa \in \{ w, a \}$ are the components, $\rho_\alpha$ is the fluid density, $S_\alpha$ is the fluid saturation, $\phi$ is the porosity of the soil, $k_{r\alpha}$ is the relative permeability for the fluid, $\mu_\alpha$ is the fluid's dynamic viscosity, $\mathbf{K}$ is the intrinsic permeability, $p_\alpha$ is the fluid pressure and $g$ is the potential of the gravity field.
+It can be derived from the two-\/phase equations, i.\-e. \[ \phi\frac{\partial S_\alpha \varrho_\alpha}{\partial t} - \text{div} \left\lbrace \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}\; \mathbf{K} \; \left( \text{\textbf{grad}} p_\alpha - \varrho_\alpha \textbf{g} \right) \right\rbrace = q_\alpha, \] where $\alpha \in \{w, n\}$ is the fluid phase, $\kappa \in \{ w, a \}$ are the components, $\rho_\alpha$ is the fluid density, $S_\alpha$ is the fluid saturation, $\phi$ is the porosity of the soil, $k_{r\alpha}$ is the relative permeability for the fluid, $\mu_\alpha$ is the fluid\textquotesingle{}s dynamic viscosity, $\mathbf{K}$ is the intrinsic permeability, $p_\alpha$ is the fluid pressure and $g$ is the potential of the gravity field.
 
 In contrast to the full two-\/phase model, the Richards model assumes gas as the non-\/wetting fluid and that it exhibits a much lower viscosity than the (liquid) wetting phase. (For example at atmospheric pressure and at room temperature, the viscosity of air is only about $1\%$ of the viscosity of liquid water.) As a consequence, the $\frac{k_{r\alpha}}{\mu_\alpha}$ term typically is much larger for the gas phase than for the wetting phase. For this reason, the Richards model assumes that $\frac{k_{rn}}{\mu_n}$ is infinitly large. This implies that the pressure of the gas phase is equivalent to the static pressure distribution and that therefore, mass conservation only needs to be considered for the wetting phase.
 
diff --git a/doc/handbook/ModelDescriptions/stokesncimplicitmodel.tex b/doc/handbook/ModelDescriptions/stokesncimplicitmodel.tex
index daa9d256500f17d4986486531a4e3e32ab2040bb..1c917d2008782326faa8afd345701bd49ce8a9d5 100644
--- a/doc/handbook/ModelDescriptions/stokesncimplicitmodel.tex
+++ b/doc/handbook/ModelDescriptions/stokesncimplicitmodel.tex
@@ -10,7 +10,7 @@ Momentum Balance\-: \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}
 
 Mass balance equation\-: \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0 \]
 
-\hyperlink{a00088}{Component} mass balance equations\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} X_g^\kappa \right) - q_g^\kappa = 0 \]
+\hyperlink{a00084}{Component} mass balance equations\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} X_g^\kappa \right) - q_g^\kappa = 0 \]
 
 This is discretized using a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit Euler method in time.
 
diff --git a/doc/handbook/ModelDescriptions/stokesncniimplicitmodel.tex b/doc/handbook/ModelDescriptions/stokesncniimplicitmodel.tex
index 410bcbfa3201ca71546cb6a860e935493867180f..a0e48a098b5f1f5da80f4c8c317f4b7b2625c570 100644
--- a/doc/handbook/ModelDescriptions/stokesncniimplicitmodel.tex
+++ b/doc/handbook/ModelDescriptions/stokesncniimplicitmodel.tex
@@ -10,7 +10,7 @@ Momentum Balance\-: \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}
 
 Mass balance equation\-: \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0 \]
 
-\hyperlink{a00088}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
+\hyperlink{a00084}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
 
 Energy balance equation\-: \[ \frac{\partial (\varrho_g u_g)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g h_g {\boldsymbol{v}}_g - \sum_\kappa \left[ h^\kappa_g D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \nabla x^\kappa_g \right] - \lambda_g \boldsymbol{\nabla} T \right) - q_T = 0 \]
 
diff --git a/doc/handbook/ModelDescriptions/zeroeqimplicitmodel.tex b/doc/handbook/ModelDescriptions/zeroeqimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3aba7e286736808c0ea850a2f994b6f14b16e4c7
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/zeroeqimplicitmodel.tex
@@ -0,0 +1,14 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements an single-\/phase isothermal free flow solving the mass and the momentum balance. For the momentum balance the Reynolds-\/averaged Navier-\/\-Stokes (R\-A\-N\-S) equation with zero equation (algebraic) turbulence model is used.
+
+Mass balance\-: \[ \frac{\partial \varrho_\textrm{g}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right) - q_\textrm{g} = 0 \]
+
+Momentum Balance\-: \[ \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}} - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right] \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g} + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right) \right) + \left(p_\textrm{g} {\bf {I}} \right) - \varrho_\textrm{g} {\bf g} = 0, \]
+
+This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
+
diff --git a/doc/handbook/ModelDescriptions/zeroeqncimplicitmodel.tex b/doc/handbook/ModelDescriptions/zeroeqncimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..d0deec2ac6103bafea1edc111d02f0b441422f66
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/zeroeqncimplicitmodel.tex
@@ -0,0 +1,16 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements an single-\/phase isothermal compositional free flow solving the mass and the momentum balance. For the momentum balance the Reynolds-\/averaged Navier-\/\-Stokes (R\-A\-N\-S) equation with zero equation (algebraic) turbulence model is used.
+
+Mass balance\-: \[ \frac{\partial \varrho_\textrm{g}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right) - q_\textrm{g} = 0 \]
+
+Momentum Balance\-: \[ \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}} - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right] \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g} + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right) \right) + \left(p_\textrm{g} {\bf {I}} \right) - \varrho_\textrm{g} {\bf g} = 0, \]
+
+\hyperlink{a00084}{Component} mass balance equations\-: \[ \frac{\partial \left(\varrho_\textrm{g} X_\textrm{g}^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} X_\textrm{g}^\kappa - \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right] \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \boldsymbol{\nabla} x_\textrm{g}^\kappa \right) - q_\textrm{g}^\kappa = 0 \]
+
+This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
+
diff --git a/doc/handbook/ModelDescriptions/zeroeqncniimplicitmodel.tex b/doc/handbook/ModelDescriptions/zeroeqncniimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c2a101598e91a6a519185b0859c679fc6eb8a467
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/zeroeqncniimplicitmodel.tex
@@ -0,0 +1,18 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements an single-\/phase non-\/isothermal compositional free flow solving the mass and the momentum balance. For the momentum balance the Reynolds-\/averaged Navier-\/\-Stokes (R\-A\-N\-S) equation with zero equation (algebraic) turbulence model is used.
+
+Mass balance\-: \[ \frac{\partial \varrho_\textrm{g}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right) - q_\textrm{g} = 0 \]
+
+Momentum Balance\-: \[ \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}} - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right] \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g} + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right) \right) + \left(p_\textrm{g} {\bf {I}} \right) - \varrho_\textrm{g} {\bf g} = 0, \]
+
+\hyperlink{a00084}{Component} mass balance equations\-: \[ \frac{\partial \left(\varrho_\textrm{g} X_\textrm{g}^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} X_\textrm{g}^\kappa - \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right] \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \boldsymbol{\nabla} x_\textrm{g}^\kappa \right) - q_\textrm{g}^\kappa = 0 \]
+
+Energy balance equation\-: \[ \frac{\partial (\varrho_\textrm{g} u_\textrm{g})}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} h_\textrm{g} {\boldsymbol{v}}_\textrm{g} - \sum_\kappa \left( h^\kappa_\textrm{g} \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right] \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \nabla x^\kappa_\textrm{g} \right) - \left[ \lambda_\textrm{g} + \lambda_\textrm{g,t} \right] \boldsymbol{\nabla} T \right) - q_\textrm{T} = 0 \]
+
+This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
+
diff --git a/doc/handbook/dumux-handbook.tex b/doc/handbook/dumux-handbook.tex
index 4d727efc9726da718a0f8874b09c8ca2aca165eb..c34e1deba7810e1aa299edf3677a3f801504b4ab 100644
--- a/doc/handbook/dumux-handbook.tex
+++ b/doc/handbook/dumux-handbook.tex
@@ -20,6 +20,7 @@
 \usepackage[hyphens]{url}
 \usepackage{tikz}
 \usepackage{pdflscape}
+\usepackage{textcomp}
 \usepackage{rotating}
 \usetikzlibrary{arrows,backgrounds,decorations.pathmorphing,patterns,positioning,fit,shapes}
 \usepackage[normalem]{ulem}
diff --git a/doc/handbook/models.tex b/doc/handbook/models.tex
index e6f1b34f3415313d8353365e2c818de5f72488d7..ce82446ed4d71ca3710c74027c10acb99ca6a1f5 100644
--- a/doc/handbook/models.tex
+++ b/doc/handbook/models.tex
@@ -249,6 +249,29 @@ subdirectories of \texttt{dumux/implicit} of the \Dumux distribution.
 \subsubsection{The Non-Isothermal $N$-Component Stokes Model: StokesNcNIModel} 
 \input{ModelDescriptions/stokesncniimplicitmodel}
 
+\subsubsection{The Isothermal Zero Equation (algebraic) Turbulence Model: ZeroeqModel} 
+\input{ModelDescriptions/zeroeqimplicitmodel}
+
+\subsubsection{The Isothermal $N$-Component Zero Equation (algebraic) Turbulence Model: ZeroeqNcModel} 
+\input{ModelDescriptions/zeroeqncimplicitmodel}
+
+\subsubsection{The Non-Isothermal $N$-Component Zero Equation (algebraic) Turbulence Model: ZeroeqNcNIModel} 
+\input{ModelDescriptions/zeroeqncniimplicitmodel}
+
+\subsubsection{The Coupled StokesNcModel - TwoPTwoCModel}
+\label{sc_2cstokes2p2c}
+\input{ModelDescriptions/2cstokes2p2cmodel}
+
+\subsubsection{The Coupled StokesNcNIModel - TwoPTwoCNIModel}
+\label{sc_2cnistokes2p2cni}
+\input{ModelDescriptions/2cnistokes2p2cnimodel}
+
+\subsubsection{The Coupled ZeroeqNcModel - TwoPTwoCModel}
+\input{ModelDescriptions/2czeroeq2p2cmodel}
+
+\subsubsection{The Coupled ZeroeqNcNIModel - TwoPTwoCNIModel}
+\input{ModelDescriptions/2cnizeroeq2p2cnimodel}
+
 \subsubsection{The Linear Elastic Model: ElasticModel} 
 \input{ModelDescriptions/elasticmodel}
 
diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index 2a3eb870a4507d53fedfc78081944fa403a72b7c..c338272db454256bb17487520bd3cfe2ef846584 100644
--- a/dumux/freeflow/stokes/stokesfluxvariables.hh
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -69,7 +69,7 @@ public:
                         const int fIdx,
                         const ElementVolumeVariables &elemVolVars,
                         const bool onBoundary = false)
-        : fvGeometry_(fvGeometry), onBoundary_(onBoundary), faceIdx_(fIdx)
+        : fvGeometry_(fvGeometry), onBoundary_(onBoundary), fIdx_(fIdx)
     {
         calculateValues_(problem, element, elemVolVars);
         determineUpwindDirection_(elemVolVars);
@@ -160,9 +160,9 @@ public:
     const SCVFace &face() const
     {
         if (onBoundary_)
-            return fvGeometry_.boundaryFace[faceIdx_];
+            return fvGeometry_.boundaryFace[fIdx_];
         else
-            return fvGeometry_.subContVolFace[faceIdx_];
+            return fvGeometry_.subContVolFace[fIdx_];
     }
 
     /*!
@@ -291,7 +291,7 @@ protected:
     // local index of the downwind vertex
     int downstreamIdx_;
     // the index of the considered face
-    int faceIdx_;
+    int fIdx_;
 };
 
 } // end namespace
diff --git a/dumux/freeflow/stokes/stokesproblem.hh b/dumux/freeflow/stokes/stokesproblem.hh
index 874a92413fc59a023491259f8f005e3e04908ed2..30a66ca70d5cd6e902807cec40a857f22c28a3a9 100644
--- a/dumux/freeflow/stokes/stokesproblem.hh
+++ b/dumux/freeflow/stokes/stokesproblem.hh
@@ -30,7 +30,8 @@
 namespace Dumux
 {
 /*!
- * \ingroup BoxStokesProblems
+ * \ingroup ImplicitBaseProblems
+ * \ingroup BoxStokesModel
  * \brief Base class for all problems which use the Stokes box model.
  *
  * This implements gravity (if desired) and a function returning the temperature.
diff --git a/dumux/freeflow/stokesnc/stokesncfluxvariables.hh b/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
index a2c99adf9b8346c6ac62d0a831e7d0ece71724a3..e0f14bb3026455ea690ff8137a17431429c20e22 100644
--- a/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
+++ b/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
@@ -36,7 +36,7 @@ namespace Dumux
 
 /*!
  * \ingroup BoxStokesncModel
- * \ingroup BoxFluxVariables
+ * \ingroup ImplicitFluxVariables
  * \brief This template class contains data which is required to
  *        calculate the component fluxes over a face of a finite
  *        volume for the compositional n component Stokes model.
diff --git a/dumux/freeflow/stokesnc/stokesncindices.hh b/dumux/freeflow/stokesnc/stokesncindices.hh
index 0c20ab0ae126de37b32823f8c6e1058797a87c24..cf21bfb2acf760dbf4a4c5b8a2fa95b425df35c9 100644
--- a/dumux/freeflow/stokesnc/stokesncindices.hh
+++ b/dumux/freeflow/stokesnc/stokesncindices.hh
@@ -34,7 +34,7 @@ namespace Dumux
 
 /*!
  * \ingroup BoxStokesncModel
- * \ingroup BoxIndices
+ * \ingroup ImplicitIndices
  * \brief The common indices for the compositional n component Stokes box model.
  *
  * \tparam PVOffset The first index in a primary variable vector.
diff --git a/dumux/freeflow/stokesnc/stokesnclocalresidual.hh b/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
index e1bd6e1457373c6f6609764260112898467d7a44..14e6e18fc7e1b1a35dd83621a5ad0c758a8492f6 100644
--- a/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
+++ b/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
@@ -35,7 +35,7 @@ namespace Dumux
 {
 /*!
  * \ingroup BoxStokesncModel
- * \ingroup BoxLocalResidual
+ * \ingroup ImplicitLocalResidual
  * \brief Element-wise calculation of the Jacobian matrix for problems
  *        using the compositional Stokes box model. This is derived
  *        from the Stokes box model.
diff --git a/dumux/freeflow/stokesnc/stokesncmodel.hh b/dumux/freeflow/stokesnc/stokesncmodel.hh
index 875bb51c3c12081f0a82aa24421fa6aebd4673ab..3c299230d1eae428ddb394d1b2451328e1417932 100644
--- a/dumux/freeflow/stokesnc/stokesncmodel.hh
+++ b/dumux/freeflow/stokesnc/stokesncmodel.hh
@@ -128,35 +128,35 @@ public:
         VolumeVariables volVars;
         ElementBoundaryTypes elemBcTypes;
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
-        for (; elemIt != endit; ++elemIt)
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
+        for (; eIt != eEndIt; ++eIt)
         {
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-            int idx = this->elementMapper().index(*elemIt);
+            int idx = this->elementMapper().index(*eIt);
 #else
-            int idx = this->elementMapper().map(*elemIt);
+            int idx = this->elementMapper().map(*eIt);
 #endif
             rank[idx] = this->gridView_().comm().rank();
 
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-            int numLocalVerts = elemIt->subEntities(dim);
+            int numLocalVerts = eIt->subEntities(dim);
 #else
-            int numLocalVerts = elemIt->template count<dim>();
+            int numLocalVerts = eIt->template count<dim>();
 #endif
             for (int i = 0; i < numLocalVerts; ++i)
             {
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-                int vIdxGlobal = this->vertexMapper().subIndex(*elemIt, i, dim);
+                int vIdxGlobal = this->vertexMapper().subIndex(*eIt, i, dim);
 #else
-                int vIdxGlobal = this->vertexMapper().map(*elemIt, i, dim);
+                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
 #endif
                 volVars.update(sol[vIdxGlobal],
                                this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                i,
                                false);
diff --git a/dumux/freeflow/stokesnc/stokesncproperties.hh b/dumux/freeflow/stokesnc/stokesncproperties.hh
index 5b15c5984f44ba0002a0c9db61f9b6fd754dd8f7..9f24e6022e88cab0180daf58711d3d3ccea5f05d 100644
--- a/dumux/freeflow/stokesnc/stokesncproperties.hh
+++ b/dumux/freeflow/stokesnc/stokesncproperties.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /*!
  * \ingroup Properties
- * \ingroup BoxProperties
+ * \ingroup ImplicitProperties
  * \ingroup BoxStokesncModel
  *
  * \file
diff --git a/dumux/freeflow/stokesnc/stokesncpropertydefaults.hh b/dumux/freeflow/stokesnc/stokesncpropertydefaults.hh
index a18dbbcd4c228b614de1f7205458c898ca4d4812..fab3f4065aa571464609b00ef6bb72249df53f35 100644
--- a/dumux/freeflow/stokesnc/stokesncpropertydefaults.hh
+++ b/dumux/freeflow/stokesnc/stokesncpropertydefaults.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /*!
  * \ingroup Properties
- * \ingroup BoxProperties
+ * \ingroup ImplicitProperties
  * \ingroup BoxStokesncModel
  *
  * \file
diff --git a/dumux/freeflow/stokesnc/stokesncvolumevariables.hh b/dumux/freeflow/stokesnc/stokesncvolumevariables.hh
index e5216e19c2e1b51fd8503b15b532a9177d090d7b..a939587f30475db761adf6aaf4d0db887e470f02 100644
--- a/dumux/freeflow/stokesnc/stokesncvolumevariables.hh
+++ b/dumux/freeflow/stokesnc/stokesncvolumevariables.hh
@@ -33,7 +33,7 @@ namespace Dumux
 
 /*!
  * \ingroup BoxStokesncModel
- * \ingroup BoxVolumeVariables
+ * \ingroup ImplicitVolumeVariables
  * \brief Contains the quantities which are are constant within a
  *        finite volume in the n-component Stokes box model.
  */
diff --git a/dumux/freeflow/stokesncni/stokesncnimodel.hh b/dumux/freeflow/stokesncni/stokesncnimodel.hh
index 857752f8a01a864d07c50302007b76f4ffe1cffd..cc30ec1b7854e36ca4cb0c73f4f0ecab5c9a48e3 100644
--- a/dumux/freeflow/stokesncni/stokesncnimodel.hh
+++ b/dumux/freeflow/stokesncni/stokesncnimodel.hh
@@ -134,35 +134,35 @@ public:
         VolumeVariables volVars;
         ElementBoundaryTypes elemBcTypes;
 		
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
-        for (; elemIt != endit; ++elemIt)
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
+        for (; eIt != eEndIt; ++eIt)
         {
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-            int idx = this->elementMapper().index(*elemIt);
+            int idx = this->elementMapper().index(*eIt);
 #else
-            int idx = this->elementMapper().map(*elemIt);
+            int idx = this->elementMapper().map(*eIt);
 #endif
             rank[idx] = this->gridView_().comm().rank();
 			
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-            int numLocalVerts = elemIt->subEntities(dim);
+            int numLocalVerts = eIt->subEntities(dim);
 #else
-            int numLocalVerts = elemIt->template count<dim>();
+            int numLocalVerts = eIt->template count<dim>();
 #endif
             for (int i = 0; i < numLocalVerts; ++i)
             {
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-                int vIdxGlobal = this->vertexMapper().subIndex(*elemIt, i, dim);
+                int vIdxGlobal = this->vertexMapper().subIndex(*eIt, i, dim);
 #else
-                int vIdxGlobal = this->vertexMapper().map(*elemIt, i, dim);
+                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
 #endif
                 volVars.update(sol[vIdxGlobal],
                                this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                i,
                                false);
diff --git a/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh b/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
index 8f94c414ec0d7cd1df72afd5383af442751aa1ac..e72aa1c50666d93c7c0c21ffe347de1fc484b947 100644
--- a/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
+++ b/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
@@ -37,8 +37,8 @@ namespace Dumux
 {
 
 /*!
- * \ingroup ImplicitFluxVariables
  * \ingroup BoxZeroEqModel
+ * \ingroup ImplicitFluxVariables
  * \brief This template class contains data which is required to
  *        calculate the component fluxes over a face of a finite
  *        volume for a ZeroEq model.
@@ -65,10 +65,10 @@ public:
     ZeroEqFluxVariables(const Problem &problem,
                         const Element &element,
                         const FVElementGeometry &fvGeometry,
-                        const int faceIdx,
+                        const int fIdx,
                         const ElementVolumeVariables &elemVolVars,
                         const bool onBoundary = false)
-        : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+        : ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
         , flowNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, FlowNormal))
         , wallNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, WallNormal))
         , eddyViscosityModel_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, EddyViscosityModel))
diff --git a/dumux/freeflow/zeroeq/zeroeqmodel.hh b/dumux/freeflow/zeroeq/zeroeqmodel.hh
index 236d205aa6560b6460d5e296e2bd39266b84e7ab..4c15939fa29c2255b8f87682ec995a3b60225203 100644
--- a/dumux/freeflow/zeroeq/zeroeqmodel.hh
+++ b/dumux/freeflow/zeroeq/zeroeqmodel.hh
@@ -35,29 +35,29 @@ namespace Dumux
  * \ingroup BoxZeroEqModel
  * \brief Adaptation of the box scheme to the ZeroEq model.
  *
- * This model implements an isothermal ZeroEq Navier Stokes flow (RANS) of a fluid
- * solving the momentum balance and the mass balance.
+ * This model implements an single-phase isothermal free flow
+ * solving the mass and the momentum balance. For the momentum balance
+ * the Reynolds-averaged Navier-Stokes (RANS) equation with zero equation
+ * (algebraic) turbulence model is used.
  *
- * \todo dynamic size of wallInterval struct
+ * Mass balance:
+ * \f[
+ *  \frac{\partial \varrho_\textrm{g}}{\partial t}
+ *  + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)
+ *  - q_\textrm{g} = 0
+ * \f]
  *
- * \todo check balance equations
  * Momentum Balance:
  * \f[
- *   \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t}
+ *   \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t}
  *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(
- *     \varrho_g {\boldsymbol{v}}_g {\boldsymbol{v}}_g}
- *   - \mu_g \left(\boldsymbol{\nabla} \boldsymbol{v}_g
- *                 + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)
- *   - \mu_{g,t} \left(\boldsymbol{\nabla} \boldsymbol{v}_g
- *                     + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)
- *   + \left(p_g {\bf {I}}
+ *     \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}}
+ *     - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right]
+ *       \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}
+ *             + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right)
  *   \right)
- *   - \varrho_g {\bf g} = 0,
- * \f]
- *
- * Mass balance:
- * \f[
- *  \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0
+ *   + \left(p_\textrm{g} {\bf {I}} \right)
+ *   - \varrho_\textrm{g} {\bf g} = 0,
  * \f]
  *
  * This is discretized by a fully-coupled vertex-centered finite volume
@@ -141,21 +141,21 @@ public:
         ElementVolumeVariables elemVolVars;
 
         // Loop over elements
-        ElementIterator elemIt = this->problem_.gridView().template begin<0>();
-        ElementIterator endit = this->problem_.gridView().template end<0>();
-        for (; elemIt != endit; ++elemIt)
+        ElementIterator eIt = this->problem_.gridView().template begin<0>();
+        ElementIterator eEndIt = this->problem_.gridView().template end<0>();
+        for (; eIt != eEndIt; ++eIt)
         {
-            if (elemIt->partitionType() != Dune::InteriorEntity)
+            if (eIt->partitionType() != Dune::InteriorEntity)
                 continue;
 
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemVolVars.update(this->problem_(), *elemIt, fvGeometry);
-            this->localResidual().evalFluxes(*elemIt, elemVolVars);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemVolVars.update(this->problem_(), *eIt, fvGeometry);
+            this->localResidual().evalFluxes(*eIt, elemVolVars);
 
             bool hasLeft = false;
             bool hasRight = false;
             for (int i = 0; i < fvGeometry.numVertices; i++) {
-                const GlobalPosition &global = fvGeometry.subContVol[i].global;
+                const GlobalPosition &globalPos = fvGeometry.subContVol[i].global;
                 if (globalI[axis] < coordVal)
                     hasLeft = true;
                 else if (globalI[axis] >= coordVal)
@@ -165,7 +165,7 @@ public:
                 continue;
 
             for (int i = 0; i < fvGeometry.numVertices; i++) {
-                const GlobalPosition &global = fvGeometry.subContVol[i].global;
+                const GlobalPosition &globalPos = fvGeometry.subContVol[i].global;
                 if (globalI[axis] < coordVal)
                     flux += this->localResidual().residual(i);
             }
@@ -213,34 +213,34 @@ public:
         {}
 
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            int idx = this->elementMapper().map(*elemIt);
+            int idx = this->elementMapper().map(*eIt);
             rank[idx] = this->gridView_().comm().rank();
 
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
 
-            int numLocalVerts = elemIt->template count<dim>();
+            int numLocalVerts = eIt->template count<dim>();
             for (int i = 0; i < numLocalVerts; ++i)
             {
-                int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
-                volVars.update(sol[globalIdx],
+                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
+                volVars.update(sol[vIdxGlobal],
                                this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                i,
                                false);
 
-                pN[globalIdx] = volVars.pressure();
-                delP[globalIdx] = volVars.pressure() - 1e5;
-                rho[globalIdx] = volVars.density();
-                mu[globalIdx] = volVars.dynamicViscosity();
-                velocity[globalIdx] = volVars.velocity();
+                pN[vIdxGlobal] = volVars.pressure();
+                delP[vIdxGlobal] = volVars.pressure() - 1e5;
+                rho[vIdxGlobal] = volVars.density();
+                mu[vIdxGlobal] = volVars.dynamicViscosity();
+                velocity[vIdxGlobal] = volVars.velocity();
             };
         }
 
@@ -254,35 +254,35 @@ public:
         // ensure that the actual values are given out
         asImp_().updateWallProperties();
 
-        elemIt = this->gridView_().template begin<0>();
-        endit = this->gridView_().template end<0>();
+        eIt = this->gridView_().template begin<0>();
+        eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            IntersectionIterator isIt = this->gridView_().ibegin(*elemIt);
-            const IntersectionIterator &endIt = this->gridView_().iend(*elemIt);
+            IntersectionIterator isIt = this->gridView_().ibegin(*eIt);
+            const IntersectionIterator &endIt = this->gridView_().iend(*eIt);
 
             for (; isIt != endIt; ++isIt)
             {
-                int faceIdx = isIt->indexInInside();
+                int fIdx = isIt->indexInInside();
 
                 FluxVariables fluxVars(this->problem_(),
-                                                    *elemIt,
+                                                    *eIt,
                                                     fvGeometry,
-                                                    faceIdx,
+                                                    fIdx,
                                                     elemVolVars,
                                                     false);
 
-                GlobalPosition globalPos = fvGeometry.subContVolFace[faceIdx].ipGlobal;
+                GlobalPosition globalPos = fvGeometry.subContVolFace[fIdx].ipGlobal;
 
                 if (asImp_().shouldWriteSCVData(globalPos))
                     asImp_().writeSCVData(volVars, fluxVars, globalPos);
@@ -760,31 +760,31 @@ public:
         FVElementGeometry fvGeometry;
         ElementBoundaryTypes elemBcTypes;
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            for (int faceIdx = 0; faceIdx < fvGeometry.numScvf; ++faceIdx)
+            for (int fIdx = 0; fIdx < fvGeometry.numScvf; ++fIdx)
             {
 
                 FluxVariables fluxVars(this->problem_(),
-                                    *elemIt,
+                                    *eIt,
                                     fvGeometry,
-                                    faceIdx,
+                                    fIdx,
                                     elemVolVars,
                                     false);
 
-                GlobalPosition globalPos = fvGeometry.subContVolFace[faceIdx].ipGlobal;
+                GlobalPosition globalPos = fvGeometry.subContVolFace[fIdx].ipGlobal;
 
                 asImp_().calculateMaxFluxVars(fluxVars, globalPos);
             }
@@ -843,34 +843,34 @@ public:
     {
         FVElementGeometry fvGeometry;
         ElementBoundaryTypes elemBcTypes;
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            IntersectionIterator isIt = this->gridView_().ibegin(*elemIt);
-            const IntersectionIterator &endIt = this->gridView_().iend(*elemIt);
+            IntersectionIterator isIt = this->gridView_().ibegin(*eIt);
+            const IntersectionIterator &endIt = this->gridView_().iend(*eIt);
 
             for (; isIt != endIt; ++isIt)
             {
-                int faceIdx = isIt->indexInInside();
+                int fIdx = isIt->indexInInside();
                 FluxVariables fluxVars(this->problem_(),
-                                                    *elemIt,
+                                                    *eIt,
                                                     fvGeometry,
-                                                    faceIdx,
+                                                    fIdx,
                                                     elemVolVars,
                                                     false);
 
-                GlobalPosition globalPos = fvGeometry.subContVolFace[faceIdx].ipGlobal;
+                GlobalPosition globalPos = fvGeometry.subContVolFace[fIdx].ipGlobal;
                 int posIdx = getPosIdx(globalPos);
                 int wallIdx = getWallIdx(globalPos, posIdx);
 
@@ -900,22 +900,22 @@ public:
         FVElementGeometry fvGeometry;
         ElementBoundaryTypes elemBcTypes;
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
+            fvGeometry.update(this->gridView_(), *eIt);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            const ReferenceElement &refElement = ReferenceElements::general(elemIt->geometry().type());
-            IntersectionIterator isIt = this->gridView_().ibegin(*elemIt);
-            const IntersectionIterator &endIt = this->gridView_().iend(*elemIt);
+            const ReferenceElement &refElement = ReferenceElements::general(eIt->geometry().type());
+            IntersectionIterator isIt = this->gridView_().ibegin(*eIt);
+            const IntersectionIterator &endIt = this->gridView_().iend(*eIt);
             for (; isIt != endIt; ++isIt)
             {
                 // handle only faces on the boundary
@@ -924,19 +924,19 @@ public:
 
                 // Assemble the boundary for all vertices of the current
                 // face
-                int faceIdx = isIt->indexInInside();
-                int numFaceVerts = refElement.size(faceIdx, 1, dim);
+                int fIdx = isIt->indexInInside();
+                int numFaceVerts = refElement.size(fIdx, 1, dim);
                 for (int faceVertIdx = 0;
                      faceVertIdx < numFaceVerts;
                      ++faceVertIdx)
                 {
-                    int boundaryFaceIdx = fvGeometry.boundaryFaceIndex(faceIdx, faceVertIdx);
+                    int boundaryFaceIdx = fvGeometry.boundaryFaceIndex(fIdx, faceVertIdx);
 
 
 
 
                     const FluxVariables boundaryVars(this->problem_(),
-                                                     *elemIt,
+                                                     *eIt,
                                                      fvGeometry,
                                                      boundaryFaceIdx,
                                                      elemVolVars,
diff --git a/dumux/freeflow/zeroeq/zeroeqproblem.hh b/dumux/freeflow/zeroeq/zeroeqproblem.hh
index ec248c572a9a5d9cca8e72a4d7881fd44061a7a0..04fa62da9947f1524c5b7da9b930cc06824e5221 100644
--- a/dumux/freeflow/zeroeq/zeroeqproblem.hh
+++ b/dumux/freeflow/zeroeq/zeroeqproblem.hh
@@ -29,7 +29,8 @@
 namespace Dumux
 {
 /*!
- * \ingroup BoxZeroEqProblems
+ * \ingroup ImplicitBaseProblems
+ * \ingroup BoxZeroEqModel
  * \brief Base class for all problems which use the ZeroEq box model.
  *
  * This implements the call of update functions used for the wall properties of
diff --git a/dumux/freeflow/zeroeqnc/zeroeqncfluxvariables.hh b/dumux/freeflow/zeroeqnc/zeroeqncfluxvariables.hh
index 84658a288e51047d95e1b78457ae05cb110a2f9d..c62f99a1b7b6b8bbd8f95b00602ac4e87db35b4c 100644
--- a/dumux/freeflow/zeroeqnc/zeroeqncfluxvariables.hh
+++ b/dumux/freeflow/zeroeqnc/zeroeqncfluxvariables.hh
@@ -37,8 +37,8 @@ namespace Dumux
 {
 
 /*!
- * \ingroup ImplicitFluxVariables
  * \ingroup BoxZeroEqncModel
+ * \ingroup ImplicitFluxVariables
  * \brief This template class contains data which is required to
  *        calculate the component fluxes over a face of a finite
  *        volume for a compositional ZeroEq model.
@@ -69,10 +69,10 @@ public:
     ZeroEqncFluxVariables(const Problem &problem,
                           const Element &element,
                           const FVElementGeometry &fvGeometry,
-                          const int faceIdx,
+                          const int fIdx,
                           const ElementVolumeVariables &elemVolVars,
                           const bool onBoundary = false)
-        : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+        : ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
         , flowNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, FlowNormal))
         , wallNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, WallNormal))
         , eddyDiffusivityModel_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, EddyDiffusivityModel))
diff --git a/dumux/freeflow/zeroeqnc/zeroeqncmodel.hh b/dumux/freeflow/zeroeqnc/zeroeqncmodel.hh
index 5ab9c1a8e5e0ae1d4076d730b454ad030f4824e3..445009d02a89c18c6690788229a900e0cdd0de29 100644
--- a/dumux/freeflow/zeroeqnc/zeroeqncmodel.hh
+++ b/dumux/freeflow/zeroeqnc/zeroeqncmodel.hh
@@ -34,33 +34,38 @@ namespace Dumux
  * \ingroup BoxZeroEqncModel
  * \brief Adaptation of the box scheme to the compositional ZeroEq model.
  *
- * This model implements an isothermal compositional ZeroEq Navier Stokes flow (RANS)
- * of a fluid solving the momentum balance,the mass balance
- * and the conservation equation for one component.
+ * This model implements an single-phase isothermal compositional free flow
+ * solving the mass and the momentum balance. For the momentum balance
+ * the Reynolds-averaged Navier-Stokes (RANS) equation with zero equation
+ * (algebraic) turbulence model is used.
  *
- * \todo update balance equations
- * Momentum Balance (has to be updated):
+ * Mass balance:
  * \f[
- *   \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t}
- *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(p_g {\bf {I}}
- *   - \mu_g \left(\boldsymbol{\nabla} \boldsymbol{v}_g
- *   + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)\right)
- *   - \varrho_g {\bf g} = 0,
+ *  \frac{\partial \varrho_\textrm{g}}{\partial t}
+ *  + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)
+ *  - q_\textrm{g} = 0
  * \f]
  *
- * \todo update balance equations
- * Mass balance (has to be updated):
+ * Momentum Balance:
  * \f[
- *  \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0
+ *   \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t}
+ *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(
+ *     \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}}
+ *     - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right]
+ *       \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}
+ *             + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right)
+ *   \right)
+ *   + \left(p_\textrm{g} {\bf {I}} \right)
+ *   - \varrho_\textrm{g} {\bf g} = 0,
  * \f]
  *
- * \todo update balance equations
- * Component mass balance equation (has to be updated):
+ * Component mass balance equations:
  * \f[
- *  \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t}
- *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa
- *   - D^\kappa_g \varrho_g \boldsymbol{\nabla} X_g^\kappa \right)
- *   - q_g^\kappa = 0
+ *  \frac{\partial \left(\varrho_\textrm{g} X_\textrm{g}^\kappa\right)}{\partial t}
+ *  + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} X_\textrm{g}^\kappa
+ *  - \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right]
+ *    \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \boldsymbol{\nabla} x_\textrm{g}^\kappa \right)
+ *  - q_\textrm{g}^\kappa = 0
  * \f]
  *
  * This is discretized by a fully-coupled vertex-centered finite volume
@@ -125,7 +130,7 @@ public:
         }
     }
 
-    //! \copydoc BoxModel::addOutputVtkFields
+    //! \copydoc ImplicitModel::addOutputVtkFields
     template <class MultiWriter>
     void addOutputVtkFields(const SolutionVector &sol,
                             MultiWriter &writer)
@@ -151,35 +156,35 @@ public:
         VolumeVariables volVars;
         ElementBoundaryTypes elemBcTypes;
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            int idx = this->elementMapper().map(*elemIt);
+            int idx = this->elementMapper().map(*eIt);
             rank[idx] = this->gridView_().comm().rank();
 
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
-            int numLocalVerts = elemIt->template count<dim>();
+            int numLocalVerts = eIt->template count<dim>();
             for (int i = 0; i < numLocalVerts; ++i)
             {
-                int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
-                volVars.update(sol[globalIdx],
+                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
+                volVars.update(sol[vIdxGlobal],
                                this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                i,
                                false);
 
-                pN[globalIdx] = volVars.pressure()*scale_;
-                delP[globalIdx] = volVars.pressure()*scale_ - 1e5;
-                Xw[globalIdx] = volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
-                rho[globalIdx] = volVars.density()*scale_*scale_*scale_;
-                mu[globalIdx] = volVars.dynamicViscosity()*scale_;
-                velocity[globalIdx] = volVars.velocity();
-                velocity[globalIdx] *= 1/scale_;
+                pN[vIdxGlobal] = volVars.pressure()*scale_;
+                delP[vIdxGlobal] = volVars.pressure()*scale_ - 1e5;
+                Xw[vIdxGlobal] = volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
+                rho[vIdxGlobal] = volVars.density()*scale_*scale_*scale_;
+                mu[vIdxGlobal] = volVars.dynamicViscosity()*scale_;
+                velocity[vIdxGlobal] = volVars.velocity();
+                velocity[vIdxGlobal] *= 1/scale_;
             }
         }
         writer.attachVertexData(pN, "P");
@@ -195,35 +200,35 @@ public:
         // just to have the actual values plotted
         asImp_().updateWallProperties();
 
-        elemIt = this->gridView_().template begin<0>();
-        endit = this->gridView_().template end<0>();
+        eIt = this->gridView_().template begin<0>();
+        eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            IntersectionIterator isIt = this->gridView_().ibegin(*elemIt);
-            const IntersectionIterator &endIt = this->gridView_().iend(*elemIt);
+            IntersectionIterator isIt = this->gridView_().ibegin(*eIt);
+            const IntersectionIterator &endIt = this->gridView_().iend(*eIt);
 
             for (; isIt != endIt; ++isIt)
             {
-                int faceIdx = isIt->indexInInside();
+                int fIdx = isIt->indexInInside();
 
                 FluxVariables fluxVars(this->problem_(),
-                                                    *elemIt,
+                                                    *eIt,
                                                     fvGeometry,
-                                                    faceIdx,
+                                                    fIdx,
                                                     elemVolVars,
                                                     false);
 
-                GlobalPosition globalPos = fvGeometry.subContVolFace[faceIdx].ipGlobal;
+                GlobalPosition globalPos = fvGeometry.subContVolFace[fIdx].ipGlobal;
 
                 if (asImp_().shouldWriteSCVData(globalPos))
                 asImp_().writeSCVData(volVars, fluxVars, globalPos);
diff --git a/dumux/freeflow/zeroeqncni/zeroeqncnifluxvariables.hh b/dumux/freeflow/zeroeqncni/zeroeqncnifluxvariables.hh
index 633391623fb8c7725a96491e39fca0170fe56300..ec6e583363cbccd4cccbf1bfb81571aa5afb7083 100644
--- a/dumux/freeflow/zeroeqncni/zeroeqncnifluxvariables.hh
+++ b/dumux/freeflow/zeroeqncni/zeroeqncnifluxvariables.hh
@@ -38,8 +38,8 @@ namespace Dumux
 {
 
 /*!
- * \ingroup ImplicitFluxVariables
  * \ingroup BoxZeroEqncniModel
+ * \ingroup ImplicitFluxVariables
  * \brief This template class contains data which is required to
  *        calculate the component fluxes over a face of a finite
  *        volume for a non-isothermal compositional ZeroEq model.
@@ -70,10 +70,10 @@ public:
     ZeroEqncniFluxVariables(const Problem &problem,
                             const Element &element,
                             const FVElementGeometry &fvGeometry,
-                            const int faceIdx,
+                            const int fIdx,
                             const ElementVolumeVariables &elemVolVars,
                             const bool onBoundary = false)
-        : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+        : ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
         , flowNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, FlowNormal))
         , wallNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, WallNormal))
         , eddyConductivityModel_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, EddyConductivityModel))
diff --git a/dumux/freeflow/zeroeqncni/zeroeqncniindices.hh b/dumux/freeflow/zeroeqncni/zeroeqncniindices.hh
index 705581158cf5188b4d4a38bf17d362ab84e1aae9..eb8e9efab6856e43a10333bd726a524a5a60e0a2 100644
--- a/dumux/freeflow/zeroeqncni/zeroeqncniindices.hh
+++ b/dumux/freeflow/zeroeqncni/zeroeqncniindices.hh
@@ -47,8 +47,8 @@ struct EddyConductivityIndices
 };
 
 /*!
- * \ingroup ImplicitIndices
  * \ingroup BoxZeroEqncniModel
+ * \ingroup ImplicitIndices
  * \brief The common indices for the non-isothermal compositional ZeroEq box model.
  *
  * \tparam PVOffset The first index in a primary variable vector.
diff --git a/dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh b/dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh
index 750541adf10f1e06224b906703952a3984af0b2f..a72abb01da382bebebf2761d73f471d8e2518781 100644
--- a/dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh
+++ b/dumux/freeflow/zeroeqncni/zeroeqncnimodel.hh
@@ -35,39 +35,48 @@ namespace Dumux
  * \ingroup BoxZeroEqncniModel
  * \brief Adaptation of the box scheme to the non-isothermal compositional ZeroEq model.
  *
- * This model implements a non-isothermal compositional ZeroEq Navier Stokes flow (RANS)
- * of a fluid solving the momentum balance,the mass balance,
- * the conservation equation for one component and the energy balance.
+ * This model implements an single-phase non-isothermal compositional free flow
+ * solving the mass and the momentum balance. For the momentum balance
+ * the Reynolds-averaged Navier-Stokes (RANS) equation with zero equation
+ * (algebraic) turbulence model is used.
  *
- * \todo update balance equations
- * Momentum Balance (has to be updated):
+ * Mass balance:
  * \f[
- *   \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t}
- *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(p_g {\bf {I}}
- *   - \mu_g \left(\boldsymbol{\nabla} \boldsymbol{v}_g
- *   + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)\right)
- *   - \varrho_g {\bf g} = 0,
+ *  \frac{\partial \varrho_\textrm{g}}{\partial t}
+ *  + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)
+ *  - q_\textrm{g} = 0
  * \f]
  *
- * \todo update balance equations
- * Mass balance (has to be updated):
+ * Momentum Balance:
  * \f[
- *  \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0
+ *   \frac{\partial \left(\varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g}\right)}{\partial t}
+ *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(
+ *     \varrho_\textrm{g} {\boldsymbol{v}_\textrm{g} {\boldsymbol{v}}_\textrm{g}}
+ *     - \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right]
+ *       \left(\boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}
+ *             + \boldsymbol{\nabla} \boldsymbol{v}_\textrm{g}^T \right)
+ *   \right)
+ *   + \left(p_\textrm{g} {\bf {I}} \right)
+ *   - \varrho_\textrm{g} {\bf g} = 0,
  * \f]
  *
- * \todo update balance equations
- * Component mass balance equation (has to be updated):
+ * Component mass balance equations:
  * \f[
- *  \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t}
- *   + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa
- *   - D^\kappa_g \varrho_g \boldsymbol{\nabla} X_g^\kappa \right)
- *   - q_g^\kappa = 0
+ *  \frac{\partial \left(\varrho_\textrm{g} X_\textrm{g}^\kappa\right)}{\partial t}
+ *  + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} X_\textrm{g}^\kappa
+ *  - \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right]
+ *    \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \boldsymbol{\nabla} x_\textrm{g}^\kappa \right)
+ *  - q_\textrm{g}^\kappa = 0
  * \f]
  *
- * \todo update balance equations
- *  Energy balance equation (has to be updated):
+ * Energy balance equation:
  * \f[
- *  -
+ *  \frac{\partial (\varrho_\textrm{g}  u_\textrm{g})}{\partial t}
+ *  + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_\textrm{g} h_\textrm{g} {\boldsymbol{v}}_\textrm{g}
+ *  - \sum_\kappa \left( h^\kappa_\textrm{g} \left[ D^\kappa_\textrm{g} + D^\kappa_\textrm{g,t} \right]
+ *                       \varrho_\textrm{g} \frac{M^\kappa}{M_\textrm{g}} \nabla x^\kappa_\textrm{g} \right)
+ *  - \left[ \lambda_\textrm{g} + \lambda_\textrm{g,t} \right] \boldsymbol{\nabla} T \right)
+ *  - q_\textrm{T} = 0
  * \f]
  *
  * This is discretized by a fully-coupled vertex-centered finite volume
@@ -132,7 +141,7 @@ public:
         }
     }
 
-    //! \copydoc BoxModel::addOutputVtkFields
+    //! \copydoc ImplicitModel::addOutputVtkFields
     template <class MultiWriter>
     void addOutputVtkFields(const SolutionVector &sol,
                             MultiWriter &writer)
@@ -159,36 +168,36 @@ public:
         VolumeVariables volVars;
         ElementBoundaryTypes elemBcTypes;
 
-        ElementIterator elemIt = this->gridView_().template begin<0>();
-        ElementIterator endit = this->gridView_().template end<0>();
+        ElementIterator eIt = this->gridView_().template begin<0>();
+        ElementIterator eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            int idx = this->elementMapper().map(*elemIt);
+            int idx = this->elementMapper().map(*eIt);
             rank[idx] = this->gridView_().comm().rank();
 
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
-            int numLocalVerts = elemIt->template count<dim>();
+            int numLocalVerts = eIt->template count<dim>();
             for (int i = 0; i < numLocalVerts; ++i)
             {
-                int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
-                volVars.update(sol[globalIdx],
+                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
+                volVars.update(sol[vIdxGlobal],
                                this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                i,
                                false);
 
-                pN[globalIdx] = volVars.pressure()*scale_;
-                delP[globalIdx] = volVars.pressure()*scale_ - 1e5;
-                Xw[globalIdx] = volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
-                rho[globalIdx] = volVars.density()*scale_*scale_*scale_;
-                mu[globalIdx] = volVars.dynamicViscosity()*scale_;
-                velocity[globalIdx] = volVars.velocity();
-                velocity[globalIdx] *= 1/scale_;
-                T[globalIdx] = volVars.temperature();
+                pN[vIdxGlobal] = volVars.pressure()*scale_;
+                delP[vIdxGlobal] = volVars.pressure()*scale_ - 1e5;
+                Xw[vIdxGlobal] = volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
+                rho[vIdxGlobal] = volVars.density()*scale_*scale_*scale_;
+                mu[vIdxGlobal] = volVars.dynamicViscosity()*scale_;
+                velocity[vIdxGlobal] = volVars.velocity();
+                velocity[vIdxGlobal] *= 1/scale_;
+                T[vIdxGlobal] = volVars.temperature();
             }
         }
         writer.attachVertexData(pN, "P");
@@ -205,35 +214,35 @@ public:
         // just to have the actual values plotted
         asImp_().updateWallProperties();
 
-        elemIt = this->gridView_().template begin<0>();
-        endit = this->gridView_().template end<0>();
+        eIt = this->gridView_().template begin<0>();
+        eEndIt = this->gridView_().template end<0>();
 
-        for (; elemIt != endit; ++elemIt)
+        for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(this->gridView_(), *elemIt);
-            elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+            fvGeometry.update(this->gridView_(), *eIt);
+            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *elemIt,
+                               *eIt,
                                fvGeometry,
                                false);
 
-            IntersectionIterator isIt = this->gridView_().ibegin(*elemIt);
-            const IntersectionIterator &endIt = this->gridView_().iend(*elemIt);
+            IntersectionIterator isIt = this->gridView_().ibegin(*eIt);
+            const IntersectionIterator &endIt = this->gridView_().iend(*eIt);
 
             for (; isIt != endIt; ++isIt)
             {
-                int faceIdx = isIt->indexInInside();
+                int fIdx = isIt->indexInInside();
 
                 FluxVariables fluxVars(this->problem_(),
-                                                    *elemIt,
+                                                    *eIt,
                                                     fvGeometry,
-                                                    faceIdx,
+                                                    fIdx,
                                                     elemVolVars,
                                                     false);
 
-                GlobalPosition globalPos = fvGeometry.subContVolFace[faceIdx].ipGlobal;
+                GlobalPosition globalPos = fvGeometry.subContVolFace[fIdx].ipGlobal;
 
                 if (asImp_().shouldWriteSCVData(globalPos))
                 asImp_().writeSCVData(volVars, fluxVars, globalPos);
diff --git a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
index 90e25e3343e3b504c325d5219d814440451c34fc..893f890f29d16b9c2602f2aa0b1ca58b3c47fe3c 100644
--- a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
+++ b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
@@ -30,8 +30,113 @@ namespace Dumux {
 
 /*!
  * \ingroup TwoPTwoCNIStokesTwoCNIModel
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  * \brief The extension of the local operator for the coupling of a two-component Stokes model
  *        and a two-phase two-component Darcy model for non-isothermal conditions.
+ *
+ * This model implements the coupling between a free-flow model
+ * and a porous-medium flow model under non-isothermal conditions.
+ * Here the coupling conditions for the individual balance are presented:
+ *
+ * The total mass balance equation:
+ * \f[
+ *  \left[
+ *    \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = -\left[
+ *      \left( \varrho_\textrm{g} \boldsymbol{v}_\textrm{g}
+ *             + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} \right) \cdot \boldsymbol{n}
+ *    \right]^\textrm{pm}
+ * \f]
+ * in which \f$n\f$ represents a vector normal to the interface pointing outside of
+ * the specified subdomain.
+ *
+ * The momentum balance (tangential), which corresponds to the Beavers-Jospeh Saffman condition:
+ * \f[
+ *  \left[
+ *   \left( {\boldsymbol{v}}_\textrm{g}
+ *          + \frac{\sqrt{\left(\boldsymbol{K} \boldsymbol{t}_i \right) \cdot \boldsymbol{t}_i}}
+ *                 {\alpha_\textrm{BJ} \mu_\textrm{g}} \boldsymbol{{\tau}}_\textrm{t} \boldsymbol{n}
+ *          \right) \cdot \boldsymbol{t}_i
+ *  \right]^\textrm{ff}
+ *  = 0
+ * \f]
+ * with
+ * \f$
+ * \boldsymbol{{\tau}_\textrm{t}} = \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right]
+ *                                  \nabla \left( \boldsymbol{v}_\textrm{g}
+ *                                                + \boldsymbol{v}_\textrm{g}^\intercal \right)
+ * \f$
+ * in which the eddy viscosity \f$ \mu_\textrm{g,t} = 0 \f$ for the Stokes equation.
+ *
+ * The momentum balance (normal):
+ * \f[
+ *  \left[
+ *    \left(
+ *      \left\lbrace
+ *        \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} {\boldsymbol{v}}_\textrm{g}^\intercal
+ *        - \boldsymbol{{\tau}}_\textrm{t}
+ *        + {p}_\textrm{g} \boldsymbol{I}
+ *      \right\rbrace \boldsymbol{n}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = p_\textrm{g}^\textrm{pm}
+ * \f]
+ *
+ * The component mass balance equation (continuity of fluxes):
+ * \f[
+ *  \left[
+ *    \left(
+ *      \varrho_\textrm{g} {X}^\kappa_\textrm{g} {\boldsymbol{v}}_\textrm{g}
+ *      - {\boldsymbol{j}}^\kappa_\textrm{g,ff,t,diff}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = -\left[
+ *    \left(
+ *      \varrho_\textrm{g} X^\kappa_\textrm{g} \boldsymbol{v}_\textrm{g}
+ *      - \boldsymbol{j}^\kappa_\textrm{g,pm,diff}
+ *      + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} X^\kappa_\textrm{l}
+ *      - \boldsymbol{j}^\kappa_\textrm{l,pm,diff}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{pm}
+ *  = 0
+ * \f]
+ * in which the diffusive fluxes \f$ j_\textrm{diff} \f$ are the diffusive fluxes as
+ * they are implemented in the individual subdomain models.
+ *
+ * The component mass balance equation (continuity of mass/ mole fractions):
+ * \f[
+ *  \left[ {X}^{\kappa}_\textrm{g} \right]^\textrm{ff}
+ *  = \left[ X^{\kappa}_\textrm{g} \right]^\textrm{pm}
+ * \f]
+ *
+ * The energy balance equation (continuity of fluxes):
+ * \f[
+ *  \left[
+ *    \left(
+ *      \varrho_\textrm{g} {h}_\textrm{g} {\boldsymbol{v}}_\textrm{g}
+ *      - {h}^\textrm{a}_\textrm{g} {\boldsymbol{j}}^\textrm{a}_\textrm{g,ff,t,diff}
+ *      - {h}^\textrm{w}_\textrm{g} {\boldsymbol{j}}^\textrm{w}_\textrm{g,ff,t,diff}
+ *      - \left( \lambda_\textrm{g} + \lambda_\textrm{g,t} \right) \nabla {T}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = -\left[
+ *       \left(
+ *         \varrho_\textrm{g} h_\textrm{g} \boldsymbol{v}_\textrm{g}
+ *         + \varrho_\textrm{l} h_\textrm{l} \boldsymbol{v}_\textrm{l}
+ *         - \lambda_\textrm{pm} \nabla T
+ *      \right) \cdot \boldsymbol{n}
+ *    \right]^\textrm{pm}
+ * \f]
+ *
+ * The energy balance equation (continuity of temperature):
+ * \f[
+ *  \left[ {T} \right]^\textrm{ff}
+ *  = \left[ T \right]^\textrm{pm}
+ * \f]
+ *
+ * This is discretized by a fully-coupled vertex-centered finite volume
+ * (box) scheme in space and by the implicit Euler method in time.
  */
 template<class TypeTag>
 class TwoCNIStokesTwoPTwoCNILocalOperator :
diff --git a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
index 3656cc2808b2c2d5ff54f7c407e06f659a38975f..6e1b70f2c8eedf42a1d86268e30a7af681072ea9 100644
--- a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
+++ b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
@@ -44,8 +44,88 @@ namespace Dumux {
 
 /*!
  * \ingroup TwoPTwoCStokesTwoCModel
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \brief The local operator for the coupling of a two-component Stokes model
  *        and a two-phase two-component porous-medium model under isothermal conditions.
+ *
+ * This model implements the coupling between a free-flow model
+ * and a porous-medium flow model under isothermal conditions.
+ * Here the coupling conditions for the individual balance are presented:
+ *
+ * The total mass balance equation:
+ * \f[
+ *  \left[
+ *    \left( \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = -\left[
+ *      \left( \varrho_\textrm{g} \boldsymbol{v}_\textrm{g}
+ *             + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} \right) \cdot \boldsymbol{n}
+ *    \right]^\textrm{pm}
+ * \f]
+ * in which \f$n\f$ represents a vector normal to the interface pointing outside of
+ * the specified subdomain.
+ *
+ * The momentum balance (tangential), which corresponds to the Beavers-Jospeh Saffman condition:
+ * \f[
+ *  \left[
+ *   \left( {\boldsymbol{v}}_\textrm{g}
+ *          + \frac{\sqrt{\left(\boldsymbol{K} \boldsymbol{t}_i \right) \cdot \boldsymbol{t}_i}}
+ *                 {\alpha_\textrm{BJ} \mu_\textrm{g}} \boldsymbol{{\tau}}_\textrm{t} \boldsymbol{n}
+ *          \right) \cdot \boldsymbol{t}_i
+ *  \right]^\textrm{ff}
+ *  = 0
+ * \f]
+ * with
+ * \f$
+ * \boldsymbol{{\tau}_\textrm{t}} = \left[ \mu_\textrm{g} + \mu_\textrm{g,t} \right]
+ *                                  \nabla \left( \boldsymbol{v}_\textrm{g}
+ *                                                + \boldsymbol{v}_\textrm{g}^\intercal \right)
+ * \f$
+ * in which the eddy viscosity \f$ \mu_\textrm{g,t} = 0 \f$ for the Stokes equation.
+ *
+ * The momentum balance (normal):
+ * \f[
+ *  \left[
+ *    \left(
+ *      \left\lbrace
+ *        \varrho_\textrm{g} {\boldsymbol{v}}_\textrm{g} {\boldsymbol{v}}_\textrm{g}^\intercal
+ *        - \boldsymbol{{\tau}}_\textrm{t}
+ *        + {p}_\textrm{g} \boldsymbol{I}
+ *      \right\rbrace \boldsymbol{n}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = p_\textrm{g}^\textrm{pm}
+ * \f]
+ *
+ * The component mass balance equation (continuity of fluxes):
+ * \f[
+ *  \left[
+ *    \left(
+ *      \varrho_\textrm{g} {X}^\kappa_\textrm{g} {\boldsymbol{v}}_\textrm{g}
+ *      - {\boldsymbol{j}}^\kappa_\textrm{g,ff,t,diff}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{ff}
+ *  = -\left[
+ *    \left(
+ *      \varrho_\textrm{g} X^\kappa_\textrm{g} \boldsymbol{v}_\textrm{g}
+ *      - \boldsymbol{j}^\kappa_\textrm{g,pm,diff}
+ *      + \varrho_\textrm{l} \boldsymbol{v}_\textrm{l} X^\kappa_\textrm{l}
+ *      - \boldsymbol{j}^\kappa_\textrm{l,pm,diff}
+ *    \right) \cdot \boldsymbol{n}
+ *  \right]^\textrm{pm}
+ *  = 0
+ * \f]
+ * in which the diffusive fluxes \f$ j_\textrm{diff} \f$ are the diffusive fluxes as
+ * they are implemented in the individual subdomain models.
+ *
+ * The component mass balance equation (continuity of mass/ mole fractions):
+ * \f[
+ *  \left[ {X}^{\kappa}_\textrm{g} \right]^\textrm{ff}
+ *  = \left[ X^{\kappa}_\textrm{g} \right]^\textrm{pm}
+ * \f]
+ *
+ * This is discretized by a fully-coupled vertex-centered finite volume
+ * (box) scheme in space and by the implicit Euler method in time.
  */
 template<class TypeTag>
 class TwoCStokesTwoPTwoCLocalOperator :
diff --git a/dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh b/dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh
index 31cdbead54b634ccb4792bf88bce996820f4f10e..f839de359713e0c04176c9031251eab21c1ceb28 100644
--- a/dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh
+++ b/dumux/multidomain/2cstokes2p2c/2cstokes2p2cnewtoncontroller.hh
@@ -32,6 +32,7 @@ namespace Dumux
 /*!
  * \ingroup Newton
  * \ingroup TwoPTwoCStokesTwoCModel
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \brief Implementation of a Newton controller for the coupling of a two-component Stokes model
  *        and a two-phase two-component porous-medium model under isothermal conditions.
  *
diff --git a/dumux/multidomain/common/multidomainproperties.hh b/dumux/multidomain/common/multidomainproperties.hh
index 22b480b770be28e83d37cdd546834b9e186aa1bd..f7a07b112ab323ef9a93afa0e61eee53c3518b50 100644
--- a/dumux/multidomain/common/multidomainproperties.hh
+++ b/dumux/multidomain/common/multidomainproperties.hh
@@ -69,12 +69,6 @@ NEW_PROP_TAG(Grid);
 //! Specifies the multidomain grid
 NEW_PROP_TAG(MultiDomainGrid);
 
-//! Specifies the time manager
-NEW_PROP_TAG(TimeManager);
-
-//! Specifies the number of equations in the coupled model
-NEW_PROP_TAG(NumEq);
-
 //! Specifies the number of equations in submodel 1
 NEW_PROP_TAG(NumEq1);
 
@@ -117,9 +111,6 @@ NEW_PROP_TAG(MultiDomainCoupling);
 NEW_PROP_TAG(MultiDomainConstraintsTrafo);
 NEW_PROP_TAG(ConstraintsTrafo);
 
-//! Specifies the type of the jacobian matrix as used for the linear solver
-NEW_PROP_TAG(JacobianMatrix);
-
 //! the routines that are used to split and merge solution vectors
 NEW_PROP_TAG(SplitAndMerge);
 
diff --git a/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
index 10ef17eaed81f79a2b4bc649e209447a9b0916f9..f11a9b18d58b1f1b00c1975885c3a786ee2b6227 100644
--- a/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
@@ -33,6 +33,7 @@ namespace Dumux
 /*!
  * \ingroup ImplicitLocalResidual
  * \ingroup TwoPTwoCStokesTwoCModel
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \brief Extending the TwoPTwoCLocalResidual by the required functions for
  *        a coupled application.
  */
diff --git a/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
index 725f1ac33d4c3e64de11afa89989f5216300374e..94bc5b549b41bd435142a2d1856be642d160a451 100644
--- a/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
@@ -35,6 +35,7 @@ namespace Dumux
 /*!
  * \ingroup ImplicitLocalResidual
  * \ingroup TwoPTwoCNIStokesTwoCNIModel
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  * \brief Extending the TwoPTwoCNILocalResidual by the required functions for
  *        a coupled application.
  */
diff --git a/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
index 52a1c32c6e69532f3601742613359e2274bb98ef..0bb98b9d7ccab4013be229e0513b0a3388aa7519 100644
--- a/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
@@ -31,6 +31,7 @@ namespace Dumux
 /*!
  * \ingroup ImplicitLocalResidual
  * \ingroup TwoPTwoCNIStokesTwoCNIModel
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  * \brief Element-wise calculation of the Jacobian matrix for problems
  *        using the coupled box model.
  */
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
index b4efd07c6677b26ae90946198e5353a6b7e6077c..6e83ffa1ea651b8a6ac2fb0aab239087c21511a3 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
@@ -33,6 +33,7 @@ namespace Dumux
 /*!
  * \ingroup ImplicitLocalResidual
  * \ingroup TwoPTwoCStokesTwoCModel
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \brief Element-wise calculation of the Jacobian matrix for problems
  *        using the coupled compositional Stokes box model.
  *        It is derived from the compositional Stokes box model.
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
index dd9c3add7b3005122d098eb40f1f9ac181f8045a..9e9d59f7de660a90d8192fff10b4e0e31315a2e8 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
@@ -34,6 +34,7 @@ namespace Dumux
 	/*!
    * \ingroup ImplicitLocalResidual
    * \ingroup TwoPTwoCNIStokesTwoCNIModel
+   * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
    * \brief Element-wise calculation of the Jacobian matrix for problems
    *        using the coupled compositional non-isothermal Stokes box model.
    *        It is derived from the compositional non-isothermal Stokes box model.
diff --git a/test/freeflow/navierstokes/navierstokestestproblem.hh b/test/freeflow/navierstokes/navierstokestestproblem.hh
index a33fac495f956af91b414f1a2e458378be38349f..0c45b1f7025741f5283f8f351250cf3b4873e333 100644
--- a/test/freeflow/navierstokes/navierstokestestproblem.hh
+++ b/test/freeflow/navierstokes/navierstokestestproblem.hh
@@ -176,8 +176,8 @@ namespace Dumux
        *        used for which equation on a given boundary segment.
        *
        * \param values The boundary types for the conservation equations
-       * \param vertex The vertex on the boundary for which the
-       *               conditions needs to be specified
+       * \param globalPos The globalPos on the boundary for which the
+       *                  conditions needs to be specified
        */
       void boundaryTypesAtPos(BoundaryTypes &values,
                             const GlobalPosition &globalPos) const
@@ -198,7 +198,7 @@ namespace Dumux
        *        control volume.
        *
        * \param values The dirichlet values for the primary variables
-       * \param vertex The vertex representing the "half volume on the boundary"
+       * \param globalPos The globalPos representing the "half volume on the boundary"
        *
        * For this method, the \a values parameter stores primary variables.
        */
diff --git a/test/freeflow/stokes/stokestestproblem.hh b/test/freeflow/stokes/stokestestproblem.hh
index 9014e768f75be22702f221b6a3e783ca814619b9..f6541a442cc2781966cd033d167ad674dfc2ffba 100644
--- a/test/freeflow/stokes/stokestestproblem.hh
+++ b/test/freeflow/stokes/stokestestproblem.hh
@@ -156,7 +156,7 @@ public:
      */
     // \{
 
-    //! \copydoc ImplicitProblem::boundaryTypes()
+    //! \copydoc ImplicitProblem::boundaryTypesAtPos()
     void boundaryTypesAtPos(BoundaryTypes &values,
                             const GlobalPosition &globalPos) const
     {
@@ -176,7 +176,7 @@ public:
             values.setDirichlet(massBalanceIdx);
     }
 
-    //! \copydoc ImplicitProblem::dirichlet()
+    //! \copydoc ImplicitProblem::dirichletAtPos()
     void dirichletAtPos(PrimaryVariables &values,
                         const GlobalPosition &globalPos) const
     {
@@ -240,7 +240,7 @@ public:
         values = Scalar(0.0);
     }
 
-    //! \copydoc ImplicitProblem::initial()
+    //! \copydoc ImplicitProblem::initialAtPos()
     void initialAtPos(PrimaryVariables &values,
                       const GlobalPosition &globalPos) const
     {
diff --git a/test/freeflow/stokes2c/stokes2ctestproblem.hh b/test/freeflow/stokes2c/stokes2ctestproblem.hh
index 0a80b032de02d6b05ca7f61951b1382bb8ef5688..f741805ad0c6148c3413d836564b6e26356c0b38 100644
--- a/test/freeflow/stokes2c/stokes2ctestproblem.hh
+++ b/test/freeflow/stokes2c/stokes2ctestproblem.hh
@@ -154,7 +154,7 @@ public:
      */
     // \{
 
-    //! \copydoc ImplicitProblem::boundaryTypes()
+    //! \copydoc ImplicitProblem::boundaryTypesAtPos()
     void boundaryTypesAtPos(BoundaryTypes &values,
                             const GlobalPosition &globalPos) const
     {
@@ -174,7 +174,7 @@ public:
             values.setDirichlet(massBalanceIdx);
     }
 
-    //! \copydoc ImplicitProblem::dirichlet()
+    //! \copydoc ImplicitProblem::dirichletAtPos()
     void dirichletAtPos(PrimaryVariables &values,
                         const GlobalPosition &globalPos) const
     {
@@ -214,7 +214,7 @@ public:
         values = Scalar(0.0);
     }
 
-    //! \copydoc ImplicitProblem::initial()
+    //! \copydoc ImplicitProblem::initialAtPos()
     void initialAtPos(PrimaryVariables &values,
                       const GlobalPosition &globalPos) const
     {
diff --git a/test/freeflow/stokes2cni/stokes2cnitestproblem.hh b/test/freeflow/stokes2cni/stokes2cnitestproblem.hh
index 3c9273351d3f8aa7ee5652ddcf0ff8d1a6cd2a5a..f8f2750d35a23be985faf82b8d11ee87b1f3b85b 100644
--- a/test/freeflow/stokes2cni/stokes2cnitestproblem.hh
+++ b/test/freeflow/stokes2cni/stokes2cnitestproblem.hh
@@ -167,7 +167,7 @@ public:
      */
     // \{
 
-    //! \copydoc ImplicitProblem::boundaryTypes()
+    //! \copydoc ImplicitProblem::boundaryTypesAtPos()
     void boundaryTypesAtPos(BoundaryTypes &values,
                             const GlobalPosition &globalPos) const
     {
@@ -188,7 +188,7 @@ public:
             values.setDirichlet(massBalanceIdx);
     }
 
-    //! \copydoc ImplicitProblem::dirichlet()
+    //! \copydoc ImplicitProblem::dirichletAtPos()
     void dirichletAtPos(PrimaryVariables &values,
                         const GlobalPosition &globalPos) const
     {
@@ -224,7 +224,7 @@ public:
         values = Scalar(0.0);
     }
 
-    //! \copydoc ImplicitProblem::initial()
+    //! \copydoc ImplicitProblem::initialAtPos()
     void initialAtPos(PrimaryVariables &values,
                       const GlobalPosition &globalPos) const
     {
diff --git a/test/freeflow/zeroeq/zeroeqchanneltestproblem.hh b/test/freeflow/zeroeq/zeroeqchanneltestproblem.hh
index 6698e21828a6dabfddc1fe81af98b8ae0932d767..30de8649836c939ea77a1250dc5a2741a5785983 100644
--- a/test/freeflow/zeroeq/zeroeqchanneltestproblem.hh
+++ b/test/freeflow/zeroeq/zeroeqchanneltestproblem.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /**
  * \file
- * \brief  Definition of an isothermal ZeroEq channel flow
+ * \brief  Definition of an isothermal ZeroEq channel flow problem.
  */
 #ifndef DUMUX_ZEROEQCHANNELTESTPROBLEM_HH
 #define DUMUX_ZEROEQCHANNELTESTPROBLEM_HH
@@ -75,11 +75,18 @@ SET_SCALAR_PROP(ZeroEqChannelTestProblem, ZeroEqWriteAllSCVData, -1);
 }
 
 /*!
+ * \ingroup BoxZeroEqModel
  * \ingroup ImplicitTestProblems
- * \ingroup ZeroEqModel
  * \brief ZeroEq problem with air flowing from the left to the right.
  *
- * This problem uses the \ref ZeroEqModel.
+ * The domain is 5.0m long and 1.5m high. Air is flow from left to right.
+ * The momentum balance has Dirichlet conditions on the left (inflow) and
+ * on bottom (no-slip), on the right it has outflow condition. On the top,
+ * which corresponds to an open surface or the center line of a pipe Neumann
+ * no-flow for tangential momentum are applied. The mass balance is outflow
+ * except at the right side.
+ *
+ * This problem uses the \ref ZeroEqModel with the modified Van Driest turbulence model.
  *
  * To run the simulation execute the following line in shell:<br>
  * <tt>./test_zeroeq_channel -ParameterFile ./test_zeroeq_channel.input</tt>
@@ -193,19 +200,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a neumann
-     *        boundary segment.
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
-     *
+     * \copydoc ImplicitProblem::neumann()
      * A neumann condition for the RANS momentum equation equation corresponds to:
-     * \f[ -\mu \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
+     * \f[ - \left[ \mu + \mu_\textrm{t} \right] \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
      */
     void neumann(PrimaryVariables &values,
                  const Element &element,
-                 const FVElementGeometry &fvElemGeom,
-                 const Intersection &is,
+                 const FVElementGeometry &fvGeometry,
+                 const Intersection &intersection,
                  int scvIdx,
                  int boundaryFaceIdx) const
     {
diff --git a/test/freeflow/zeroeq/zeroeqtestproblem.hh b/test/freeflow/zeroeq/zeroeqtestproblem.hh
index 087917a9d14ccfe6bd48b3cd93a8f5d97e53839e..a85845fac78830b190efaff9da214b5a826ad62a 100644
--- a/test/freeflow/zeroeq/zeroeqtestproblem.hh
+++ b/test/freeflow/zeroeq/zeroeqtestproblem.hh
@@ -18,10 +18,10 @@
  *****************************************************************************/
 /**
  * \file
- * \brief  Definition of an isothermal ZeroEq problem.
+ * \brief Definition of an isothermal ZeroEq problem.
  *
- * This problem implements an experiment performed by John Laufer
- * (The structure of turbulence in fully developed pipe flow, NACA Report, 1954).
+ * This problem implements a pipe flow experiment performed by John Laufer
+ * (Laufer J., The structure of turbulence in fully developed pipe flow, NACA Report, 1954).
  */
 #ifndef DUMUX_ZEROEQTESTPROBLEM_HH
 #define DUMUX_ZEROEQTESTPROBLEM_HH
@@ -75,19 +75,22 @@ SET_SCALAR_PROP(ZeroEqTestProblem, ZeroEqWriteAllSCVData, .8875);
 }
 
 /*!
+ * \ingroup BoxZeroEqModel
  * \ingroup ImplicitTestProblems
- * \ingroup ZeroEqModel
  * \brief ZeroEq problem with air flowing from the left to the right.
  *
  * The domain is sized 10m times 0.2469m. The problem is taken from an
- * experimental setup by John Laufer (The structure of turbulence in fully developed pipe flow,
- * NACA Report, 1954). The boundary conditions for the momentum balances
+ * experimental setup by John Laufer (J. Laufer, The structure of turbulence in
+ * fully developed pipe flow, NACA Report, 1954).
+ * The boundary conditions for the momentum balances
  * are set to Dirichlet on the left (inflow) and outflow on the right boundary.
  * The mass balance has outflow boundary conditions, which are replaced in the
- * localresidual by the sum of the two momentum balances. On the right boundary
+ * localresidual by the sum of the two momentum balances. On the right boundary,
  * the mass balance receives a Dirichlet value to set the pressure level.
  *
- * This problem uses the \ref ZeroEqModel.
+ * This problem uses the \ref ZeroEqModel with the Baldwin-Lomax turbulence model
+ * (Baldwin, B. S. & Lomax, H. Thin Layer Approximation and Algebraic Model for
+ * Seperated Turbulent Flows AIAA Journal, 1978).
  *
  * To run the simulation execute the following line in shell:<br>
  * <tt>./test_zeroeq -ParameterFile ./test_zeroeq.input</tt>
@@ -193,19 +196,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a neumann
-     *        boundary segment.
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
-     *
+     * \copydoc ImplicitProblem::neumann()
      * A neumann condition for the RANS momentum equation equation corresponds to:
-     * \f[ -\mu \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
+     * \f[ - \left[ \mu + \mu_\textrm{t} \right] \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
      */
     void neumann(PrimaryVariables &values,
                  const Element &element,
-                 const FVElementGeometry &fvElemGeom,
-                 const Intersection &is,
+                 const FVElementGeometry &fvGeometry,
+                 const Intersection &intersection,
                  int scvIdx,
                  int boundaryFaceIdx) const
     {
diff --git a/test/freeflow/zeroeq2c/zeroeq2ctestproblem.hh b/test/freeflow/zeroeq2c/zeroeq2ctestproblem.hh
index 70eb9884d5134d91676982be2cd36768aebd4497..226d58a1b09c37f000f29840ec63ef866de7d13d 100644
--- a/test/freeflow/zeroeq2c/zeroeq2ctestproblem.hh
+++ b/test/freeflow/zeroeq2c/zeroeq2ctestproblem.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /**
  * \file
- * \brief  Definition of an isothermal compositional ZeroEq problem
+ * \brief  Definition of an isothermal compositional ZeroEqnc problem
  */
 #ifndef DUMUX_ZEROEQTWOCTESTPROBLEM_HH
 #define DUMUX_ZEROEQTWOCTESTPROBLEM_HH
@@ -69,14 +69,19 @@ SET_BOOL_PROP(ZeroEq2cTestProblem, ProblemEnableGravity, false);
 }
 
 /*!
+ * \ingroup BoxZeroEqnCModel
  * \ingroup ImplicitTestProblems
- * \ingroup ZeroEqTwoCModel
- * \brief Compositional ZeroEq flow problem .
+ * \brief isothermal compositional ZeroEqnc flow problem with injection from the left.
  *
- * \todo This is just some arbitrary test problem, for which the description
- *       has to be added
+ * The domain is sized 3m times 1m. Air is flowing in a pipe from left to right.
+ * Thus the momentum balance has Dirichlet conditions (inflow) on the left and
+ * no-slip on top and on bottom, on the right outflow conditions are applied. The total
+ * mass balance has outflow boundary conditions everywhere, except at the right where
+ * Dirichlet values are set. Water vapor is injected in the middle of the left border
+ * and is mixed and transported in the air stream.
  *
- * This problem uses the \ref ZeroEqTwoCModel.
+ * This problem uses the \ref ZeroEqncModel with the Prandtl mixing length model for
+ * the eddy viscosity and the model by Meier and Rotta for the eddy diffusivity.
  *
  * To run the simulation execute the following line in shell:<br>
  * <tt>./test_zeroeq2c -ParameterFile ./test_zeroeq2c.input</tt>
@@ -157,17 +162,14 @@ public:
     // \{
 
     /*!
-     * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given boundary segment.
-     *
-     * \param values The boundary types for the conservation equations
-     * \param vertex The vertex on the boundary for which the
-     *               conditions needs to be specified
+     * \name Boundary conditions
      */
-    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
-    {
-        const GlobalPosition globalPos = vertex.geometry().center();
+    // \{
 
+    //! \copydoc ImplicitProblem::boundaryTypesAtPos()
+    void boundaryTypesAtPos(BoundaryTypes &values,
+                            const GlobalPosition &globalPos) const
+    {
         values.setAllDirichlet();
 
         if (onRightBoundary_(globalPos)
@@ -190,19 +192,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a neumann
-     *        boundary segment.
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
-     *
-     * A NEUMANN condition for the Stokes equation corresponds to:
-     * \f[ -\mu \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
+     * \copydoc ImplicitProblem::neumann()
+     * A neumann condition for the RANS momentum equation equation corresponds to:
+     * \f[ - \left[ \mu + \mu_\textrm{t} \right] \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
      */
     void neumann(PrimaryVariables &values,
                  const Element &element,
-                 const FVElementGeometry &fvElemGeom,
-                 const Intersection &is,
+                 const FVElementGeometry &fvGeometry,
+                 const Intersection &intersection,
                  int scvIdx,
                  int boundaryFaceIdx) const
     {
diff --git a/test/freeflow/zeroeq2cni/zeroeq2cnitestproblem.hh b/test/freeflow/zeroeq2cni/zeroeq2cnitestproblem.hh
index 85c66e36e4adb4e117dd3a69b2352673730fcd94..678f782b11061e31cf8ba0470fd625fc4dcb2d05 100644
--- a/test/freeflow/zeroeq2cni/zeroeq2cnitestproblem.hh
+++ b/test/freeflow/zeroeq2cni/zeroeq2cnitestproblem.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /**
  * \file
- * \brief  Definition of a simple non-isothermal compositional ZeroEq2cni box problem
+ * \brief Definition of a non-isothermal compositional ZeroEqncni problem
  */
 #ifndef DUMUX_ZEROEQTWOCNITESTPROBLEM_HH
 #define DUMUX_ZEROEQTWOCNITESTPROBLEM_HH
@@ -69,14 +69,19 @@ SET_BOOL_PROP(ZeroEq2cniTestProblem, ProblemEnableGravity, false);
 }
 
 /*!
- * \ingroup ZeroEq2cniModel
+ * \ingroup BoxZeroEqncniModel
  * \ingroup ImplicitTestProblems
- * \brief Non-isothermal compositional ZeroEq2cni flow problem.
+ * \brief Non-isothermal compositional ZeroEqncni flow problem.
  *
- * \todo This is just some arbitrary test problem, for which the description
- *       has to be added
- * 
- * This problem uses the \ref BoxZeroEq2cniModel.
+ * The domain is sized 3m times 1m. Air is flowing in a pipe from left to right.
+ * Thus the momentum balance has Dirichlet conditions (inflow) on the left and
+ * no-slip on top and on bottom, on the right outflow conditions are applied. The total
+ * mass balance has outflow boundary conditions everywhere, except at the right where
+ * Dirichlet values are set. The energy equation has Dirichlet conditions everywhere,
+ * except on the right. The bottom of the pipe is cooler than the fluid.
+ *
+ * This problem uses the \ref ZeroEqncniModel with the Prandtl mixing length model for
+ * the eddy viscosity and the model by Deissler for the eddy diffusivity and eddy conductivity.
  *
  * To run the simulation execute the following line in shell:
  * <tt>./test_zeroeq2cni -ParameterFile ./test_zeroeq2cni.input</tt>
@@ -179,19 +184,14 @@ public:
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a neumann
-     *        boundary segment.
-     *
-     * For this method, the \a values parameter stores the mass flux
-     * in normal direction of each phase. Negative values mean influx.
-     *
-     * A NEUMANN condition for the Stokes equation corresponds to:
-     * \f[ -\mu \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
+     * \copydoc ImplicitProblem::neumann()
+     * A neumann condition for the RANS momentum equation equation corresponds to:
+     * \f[ - \left[ \mu + \mu_\textrm{t} \right] \nabla {\bf v} \cdot {\bf n} + p \cdot {\bf n} = q_N \f]
      */
     void neumann(PrimaryVariables &values,
                  const Element &element,
-                 const FVElementGeometry &fvElemGeom,
-                 const Intersection &is,
+                 const FVElementGeometry &fvGeometry,
+                 const Intersection &intersection,
                  int scvIdx,
                  int boundaryFaceIdx) const
     {
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
index 335e3509e0340a0a22781d34be1b58364f565535..088cf801957b41fad73c8545069943e65ffd33c0 100644
--- a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
@@ -110,7 +110,7 @@ SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNIProblem, LinearSolver, SuperLUBackend<TypeTa
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCNIStokesTwoCNIModel
  * \brief The problem class for the coupling of a non-isothermal two-component Stokes
  *        and a non-isothermal two-phase two-component Darcy model.
  *
@@ -294,22 +294,22 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        ElementIterator eendit = mdGrid.template leafend<0>();
-        for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
-             elementIt != eendit; ++elementIt)
+        ElementIterator eEndIt = mdGrid.template leafend<0>();
+        for (ElementIterator eIt = mdGrid.template leafbegin<0>();
+             eIt != eEndIt; ++eIt)
         {
             // this is required for parallelization
             // checks if element is within a partition
-            if (elementIt->partitionType() != Dune::InteriorEntity)
+            if (eIt->partitionType() != Dune::InteriorEntity)
                 continue;
 
-            GlobalPosition globalPos = elementIt->geometry().center();
+            GlobalPosition globalPos = eIt->geometry().center();
 
             if (globalPos[1] > interfacePosY_)
-                mdGrid.addToSubDomain(stokes2cni_,*elementIt);
+                mdGrid.addToSubDomain(stokes2cni_,*eIt);
             else
                 if(globalPos[0] > noDarcyX_)
-                    mdGrid.addToSubDomain(twoPtwoCNI_,*elementIt);
+                    mdGrid.addToSubDomain(twoPtwoCNI_,*eIt);
         }
         mdGrid.preUpdateSubDomains();
         mdGrid.updateSubDomains();
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
index 63af47ccabcf5e7b85f9637c3a52b5f1a43f4780..393cc7823f0a626e159b17c437761beb921fc871 100644
--- a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
@@ -59,7 +59,7 @@ SET_TYPE_PROP(TwoCNIStokesTwoPTwoCNISpatialParams,
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCNIStokesTwoCNIModel
  * \brief Definition of the spatial parameters for
  *        the coupling of an non-isothermal two-component Stokes
  *        and an non-isothermal two-phase two-component Darcy model.
diff --git a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
index 3d84c0e69e9dc521b5e39654fd8386bb54ffc3bd..048e06dc8cab0eb62f6307aad123dc3e55884978 100644
--- a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
@@ -83,7 +83,7 @@ SET_BOOL_PROP(TwoPTwoCNISubProblem, ProblemEnableGravity, true);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCNIStokesTwoCNIModel
  * \brief Non-isothermal two-phase two-component porous-medium subproblem
  *        with coupling at the top boundary.
  *
@@ -202,7 +202,7 @@ public:
                 << "WaterMass"
                 << std::endl;
     }
-
+   //! \brief The destructor
     ~TwoPTwoCNISubProblem()
     {
         evaporationFile.close();
@@ -345,12 +345,12 @@ public:
     /*!
      * \brief Return the initial phase state inside a control volume.
      *
-     * \param vert The vertex
-     * \param dofIdxGlobal The global index of the degree of freedom
+     * \param vertex The vertex
+     * \param vIdxGlobal The global index of the vertex
      * \param globalPos The global position
      */
-    int initialPhasePresence(const Vertex &vert,
-                             const int &dofIdxGlobal,
+    int initialPhasePresence(const Vertex &vertex,
+                             const int &vIdxGlobal,
                              const GlobalPosition &globalPos) const
     {
         return bothPhases;
diff --git a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
index 319a1b894283bd56054355604ef7ec7ea66d2812..2e72667623119357092ba4d9d049924eb57e3918 100644
--- a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
@@ -72,7 +72,7 @@ SET_BOOL_PROP(Stokes2cniSubProblem, EnableNavierStokes, false);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCNIStokesTwoCNIModel
  * \brief Non-isothermal two-component stokes subproblem with air flowing
  *        from the left to the right and coupling at the bottom.
  *
diff --git a/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
index c27d8754ce33ecddf5fa7ddb772efa484b789384..64abdff3a38497d60057c66028a77b957622a47c 100644
--- a/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
+++ b/test/multidomain/2cnistokes2p2cni/test_2cnistokes2p2cni.cc
@@ -130,7 +130,7 @@ int startLocal_(int argc, char **argv,
 
     std::string dgfFileName;
     Scalar dt, tEnd;
-    Dune::FieldVector<int, dim> nElements;
+    Dune::FieldVector<int, dim> numElements;
     Scalar interfacePosY, gradingFactorY;
     int gridRefinement;
     bool useInterfaceMeshCreator;
@@ -138,8 +138,8 @@ int startLocal_(int argc, char **argv,
     dgfFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Grid, File);
     dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
     tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
-	nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
-    nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+    numElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+    numElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
     interfacePosY = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
     gradingFactorY = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingFactorY);
     gridRefinement = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, Refinement);
@@ -150,7 +150,7 @@ int startLocal_(int argc, char **argv,
     if (useInterfaceMeshCreator)
     {
         Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
-        GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePosY, gradingFactorY);
+        GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, numElements, interfacePosY, gradingFactorY);
     }
     else
     {
@@ -240,6 +240,7 @@ int startLocal_(int argc, char **argv,
  *
  * \param argc  The number of command line arguments of the program
  * \param argv  The contents of the command line arguments of the program
+ * \param printUsage  Print a usage string for simulations.
  */
 template <class TypeTag>
 int startLocal(int argc, char **argv,
diff --git a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
index ef4ed40017ca5904f0a6e6fc1017d8018d5c1a1f..821592ae53f03ca04e185377f39c14e1cf73f427 100644
--- a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
+++ b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /**
  * \file
- * \brief The problem class for the coupling of a non-isothermal two-component ZeroEq
+ * \brief The problem which couples a non-isothermal two-component ZeroEq
  *        and a non-isothermal two-phase two-component Darcy model.
  */
 #ifndef DUMUX_TWOCNIZEROEQTWOPTWOCNIPROBLEM_HH
@@ -96,15 +96,15 @@ SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNIProblem, LinearSolver, UMFPackBackend<TypeTa
 }
 
 /*!
- * \brief The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2cni)
- *        and a non-isothermal two-phase two-component Darcy model (2p2cni).
+ * \ingroup ImplicitTestProblems
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  *
- *        The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2cni)
+ * \brief The problem which couples a non-isothermal two-component ZeroEq (zeroeq2cni)
  *        and a non-isothermal two-phase two-component Darcy model (2p2cni).
- *        It uses the 2p2cniCoupling model and the ZeroEq2cnicoupling model and provides
- *        the problem specifications for common parameters of the two submodels.
- *        The initial and boundary conditions of the submodels are specified in the two subproblems,
- *        2p2cnisubproblem.hh and zeroeq2cnisubproblem.hh, which are accessible via the coupled problem.
+ *
+ * It uses the multidomain problem and specifies parameters for the two submodels.
+ * The initial and boundary conditions of the submodels are specified in the two subproblems,
+ * 2p2csubproblem.hh and zeroeq2csubproblem.hh, which are accessible via the coupled problem.
  */
 template <class TypeTag = TTAG(TwoCNIZeroEqTwoPTwoCNIProblem) >
 class TwoCNIZeroEqTwoPTwoCNIProblem : public MultiDomainProblem<TypeTag>
@@ -125,9 +125,9 @@ class TwoCNIZeroEqTwoPTwoCNIProblem : public MultiDomainProblem<TypeTag>
 
 public:
     /*!
-     * \brief docme
+     * \brief The problem for the coupling of ZeroEq and Darcy flow
      *
-     * \param hostGrid docme
+     * \param mdGrid The multidomain grid
      * \param timeManager The TimeManager which is used by the simulation
      *
      */
@@ -172,21 +172,21 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        ElementIterator eendit = mdGrid.template leafend<0>();
-        for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
-             elementIt != eendit; ++elementIt)
+        ElementIterator eEndIt = mdGrid.template leafend<0>();
+        for (ElementIterator eIt = mdGrid.template leafbegin<0>();
+             eIt != eEndIt; ++eIt)
         {
             // this is required for parallelization checks if element is within a partition
-            if (elementIt->partitionType() != Dune::InteriorEntity)
+            if (eIt->partitionType() != Dune::InteriorEntity)
                 continue;
 
-            GlobalPosition globalPos = elementIt->geometry().center();
+            GlobalPosition globalPos = eIt->geometry().center();
 
             if (globalPos[1] > interfacePos_)
-                mdGrid.addToSubDomain(zeroeq2cni_,*elementIt);
+                mdGrid.addToSubDomain(zeroeq2cni_,*eIt);
             else
                 if(globalPos[0] > noDarcyX1_ && globalPos[0] < noDarcyX2_)
-                    mdGrid.addToSubDomain(twoPtwoCNI_,*elementIt);
+                    mdGrid.addToSubDomain(twoPtwoCNI_,*eIt);
         }
         mdGrid.preUpdateSubDomains();
         mdGrid.updateSubDomains();
@@ -196,11 +196,11 @@ public:
         gridinfo(this->sdGrid2());
     }
 
-    //! \copydoc Dumux::CoupledProblem::episodeEnd()
+    //! \copydoc Dumux::ImplicitProblem::episodeEnd()
     void episodeEnd()
     { this->timeManager().startNextEpisode(episodeLength_); }
 
-    //! \copydoc Dumux::CoupledProblem::shouldWriteRestartFile()
+    //! \copydoc Dumux::ImplicitProblem::shouldWriteRestartFile()
     bool shouldWriteRestartFile() const
     {
         return ( ((this->timeManager().timeStepIndex() > 0)
@@ -209,7 +209,7 @@ public:
                 || this->timeManager().episodeWillBeOver());
     }
 
-    //! \copydoc Dumux::CoupledProblem::shouldWriteOutput()
+    //! \copydoc Dumux::ImplicitProblem::shouldWriteOutput()
     bool shouldWriteOutput() const
     {
         return ( ((this->timeManager().timeStepIndex() > 0)
diff --git a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh
index 20e8c5e408dc3403000c2317432a11ea27a0bdec..daa8aaaecf03899ddd0b38debdd7db54f858601d 100644
--- a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh
+++ b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cnispatialparameters.hh
@@ -58,9 +58,8 @@ SET_TYPE_PROP(TwoCNIZeroEqTwoPTwoCNISpatialParams,
 
 
 /*!
- * \ingroup TwoPTwoCNIModel
- * \ingroup ZeroEqTwoCNIModel
  * \ingroup ImplicitTestProblems
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  * \brief Definition of the spatial parameters for
  *        the coupling of a non-isothermal two-component ZeroEq
  *        and a non-isothermal two-phase two-component Darcy model.
diff --git a/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh
index 8412a9258097c0f51b85e4c007c2d88b8d0e7978..0098565bdc4f3862a46e02157ec49ae9f81ede11 100644
--- a/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh
+++ b/test/multidomain/2cnizeroeq2p2cni/2p2cnisubproblem.hh
@@ -18,7 +18,6 @@
  *****************************************************************************/
 /*!
  * \file
- *
  * \brief Non-isothermal two-phase two-component porous-medium subproblem
  *        with coupling at the top boundary.
  */
@@ -82,14 +81,18 @@ SET_BOOL_PROP(TwoPTwoCNISubProblem, ProblemEnableGravity, true);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
  * \brief Non-isothermal two-phase two-component porous-medium subproblem
  *        with coupling at the top boundary.
  *
- * \todo update description
+ * The porous-medium subdomain is sized 0.25m times 0.25m. The boundary conditions
+ * are Neumann no-flow everywhere, except at the top, where coupling conditions
+ * are applied to all balance equations. They handle the exchange to the free-flow
+ * subdomain. At the bottom of the porous-medium subdomain a constant temperature is
+ * set.
  *
- * This sub problem uses the \ref TwoPTwoCModel. It is part of the 2p2cni model and
- * is combined with the zeroeq2cnisubproblem for the free flow domain.
+ * This subproblem uses the \ref TwoPTwoCNIModel. It is part of a multidomain model and
+ * combined with the zeroeq2cnisubproblem for the free flow domain.
  */
 template <class TypeTag = TTAG(TwoPTwoCNISubProblem) >
 class TwoPTwoCNISubProblem : public ImplicitPorousMediaProblem<TypeTag>
@@ -275,11 +278,11 @@ public:
      * \brief Return the initial phase state inside a control volume.
      *
      * \param vertex The vertex
-     * \param globalIdx The index of the global vertex
+     * \param vIdxGlobal The global index of the vertex
      * \param globalPos The global position
      */
     int initialPhasePresence(const Vertex &vertex,
-                             const int &globalIdx,
+                             const int &vIdxGlobal,
                              const GlobalPosition &globalPos) const
     {
         return bothPhases;
diff --git a/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc
index 82557555a92df481132a441028a3abeddfa5cdce..2d55ac65a08270387b66d0d46beb6cc1f2a7d4f9 100644
--- a/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc
+++ b/test/multidomain/2cnizeroeq2p2cni/test_2cnizeroeq2p2cni.cc
@@ -173,7 +173,7 @@ int startLocal_(int argc, char **argv,
 
     std::string dgfFileName = "temp.dgf";
     Scalar dt, tEnd;
-    Dune::FieldVector<int, dim> nElements;
+    Dune::FieldVector<int, dim> numElements;
     Scalar interfacePos, gradingFactor;
     bool refineTop;
 
@@ -182,8 +182,8 @@ int startLocal_(int argc, char **argv,
         dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
         tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
 
-        nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
-        nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+        numElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+        numElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
         interfacePos = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
         gradingFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingY);
         refineTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, RefineTop);
@@ -237,7 +237,7 @@ int startLocal_(int argc, char **argv,
     }
 
     Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
-    GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePos, gradingFactor, refineTop);
+    GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, numElements, interfacePos, gradingFactor, refineTop);
 
     if (mpiHelper.size() > 1) {
         if (!Dune::Capabilities::isParallel<Grid>::v) {
@@ -311,6 +311,7 @@ int startLocal_(int argc, char **argv,
  *
  * \param argc  The number of command line arguments of the program
  * \param argv  The contents of the command line arguments of the program
+ * \param printUsage  Print a usage string for simulations.
  */
 template <class TypeTag>
 int startLocal(int argc, char **argv,
diff --git a/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh b/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh
index c045774b3ce16c2169a992b25d07aed9da0c8536..2859288373a0e309cb1d1e91146783e6fd184b46 100644
--- a/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh
+++ b/test/multidomain/2cnizeroeq2p2cni/zeroeq2cnisubproblem.hh
@@ -87,12 +87,17 @@ SET_SCALAR_PROP(ZeroEq2cniSubProblem, FreeFlowSinusTemperaturePeriod, 3600.0);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
- * \brief ZeroEq2cni problem with air flowing from the left to the right.
+ * \ingroup TwoPTwoCNIZeroEqTwoCNIModel
+ * \brief Non-isothermal two-component ZeroEq subproblem with air flowing
+ *        from the left to the right and coupling at the bottom.
+ *
+ * The free-flow subdomain is sized 0.75m times 0.5m. Dry and hot air is flowing from left (Dirichlet)
+ * to right (outflow), at the middle third of the bottom the coupling conditions
+ * are applied to all balance equations. They handle the exchange to the porous-medium
+ * subdomain.
  *
- * \todo update test description
- * This sub problem uses the \ref ZeroEq2cniModel. It is part of the 2cnizeroeq2p2cni model and
- * is combined with the 2p2csubproblem for the Darcy domain.
+ * This subproblem uses the \ref ZeroEqncniModel. It is part of a multidomain model and
+ * combined with the 2p2cnisubproblem for the porous-medium domain.
  */
 template <class TypeTag>
 class ZeroEq2cniSubProblem : public ZeroEqProblem<TypeTag>
@@ -173,11 +178,11 @@ public:
     }
 
     // functions have to be overwritten, otherwise they remain uninitialised
-    //! \copydoc BoxProblem::&bBoxMin()
+    //! \copydoc Dumux::ImplicitProblem::bBoxMin()
     const GlobalPosition &bBoxMin() const
     { return bBoxMin_; }
 
-    //! \copydoc BoxProblem::&bBoxMax()
+    //! \copydoc Dumux::ImplicitProblem::bBoxMax()
     const GlobalPosition &bBoxMax() const
     { return bBoxMax_; }
 
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
index 6861b97ee2e4ce75b4134482a4f3e879544cf4fe..721a6c9df2dcd12b23eea4c5d0ef282c1f082758 100644
--- a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
@@ -110,7 +110,7 @@ SET_TYPE_PROP(TwoCStokesTwoPTwoCProblem, LinearSolver, SuperLUBackend<TypeTag>);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCStokesTwoCModel
  * \brief The problem class for the coupling of an isothermal two-component Stokes
  *        and an isothermal two-phase two-component Darcy model.
  *
@@ -285,22 +285,22 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        ElementIterator eendit = mdGrid.template leafend<0>();
-        for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
-             elementIt != eendit; ++elementIt)
+        ElementIterator eEndit = mdGrid.template leafend<0>();
+        for (ElementIterator eIt = mdGrid.template leafbegin<0>();
+             eIt != eEndit; ++eIt)
         {
             // this is required for parallelization
             // checks if element is within a partition
-            if (elementIt->partitionType() != Dune::InteriorEntity)
+            if (eIt->partitionType() != Dune::InteriorEntity)
                 continue;
 
-            GlobalPosition globalPos = elementIt->geometry().center();
+            GlobalPosition globalPos = eIt->geometry().center();
 
             if (globalPos[1] > interfacePosY_)
-                mdGrid.addToSubDomain(stokes2c_,*elementIt);
+                mdGrid.addToSubDomain(stokes2c_,*eIt);
             else
                 if(globalPos[0] > noDarcyX_)
-                    mdGrid.addToSubDomain(twoPtwoC_,*elementIt);
+                    mdGrid.addToSubDomain(twoPtwoC_,*eIt);
         }
         mdGrid.preUpdateSubDomains();
         mdGrid.updateSubDomains();
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
index 052f371c2baed3edbd872f3cc5b7371a84ed1b23..a72aa09b4b17f8c989d6b864850aa095d301a208 100644
--- a/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
@@ -60,7 +60,7 @@ SET_TYPE_PROP(TwoCStokesTwoPTwoCSpatialParams,
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCStokesTwoCModel
  * \brief Definition of the spatial parameters for
  *        the coupling of an isothermal two-component Stokes
  *        and an isothermal two-phase two-component Darcy model.
diff --git a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
index ca8b208150af0ed6905b8c9d8f14f75e6ef54b74..9d7d6b9de1a6b8398f92ad35f97002c78052baa7 100644
--- a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
@@ -77,7 +77,7 @@ SET_BOOL_PROP(TwoPTwoCSubProblem, ProblemEnableGravity, true);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCStokesTwoCModel
  * \brief Isothermal two-phase two-component porous-medium subproblem
  *        with coupling at the top boundary.
  *
@@ -318,11 +318,11 @@ public:
      * \brief Return the initial phase state inside a control volume.
      *
      * \param vertex The vertex
-     * \param dofIdxGlobal The global index of the degree of freedom
+     * \param vIdxGlobal The global index of the vertex
      * \param globalPos The global position
      */
     int initialPhasePresence(const Vertex &vertex,
-                             const int &dofIdxGlobal,
+                             const int &vIdxGlobal,
                              const GlobalPosition &globalPos) const
     {
         return bothPhases;
diff --git a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
index 41232124c9fb68accefcab61dbd95fe82e036312..1af3c08b0c05c4b3ff67371dfb07dbe07b295c0e 100644
--- a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
@@ -69,7 +69,7 @@ SET_BOOL_PROP(Stokes2cSubProblem, EnableNavierStokes, false);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCStokesTwoCModel
  * \brief Isothermal two-component stokes subproblem with air flowing
  *        from the left to the right and coupling at the bottom.
  *
diff --git a/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
index 40c3681b6b85be23360b3e6aa820124076ac4b61..7d6483d9993389bbb3b5d5ae94f75fe82bd8a8e4 100644
--- a/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
+++ b/test/multidomain/2cstokes2p2c/test_2cstokes2p2c.cc
@@ -130,15 +130,15 @@ int startLocal_(int argc, char **argv,
 
     std::string dgfFileName;
     Scalar dt, tEnd;
-    Dune::FieldVector<int, dim> nElements;
+    Dune::FieldVector<int, dim> numElements;
     Scalar interfacePosY, gradingFactorY;
     bool useInterfaceMeshCreator;
 
     dgfFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Grid, File);
     dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
     tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
-    nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
-    nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+    numElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+    numElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
     interfacePosY = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePosY);
     gradingFactorY = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingFactorY);
     useInterfaceMeshCreator = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, UseInterfaceMeshCreator);
@@ -148,7 +148,7 @@ int startLocal_(int argc, char **argv,
     if (useInterfaceMeshCreator)
     {
         Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
-        GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePosY, gradingFactorY);
+        GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, numElements, interfacePosY, gradingFactorY);
     }
     else
     {
@@ -235,6 +235,7 @@ int startLocal_(int argc, char **argv,
  *
  * \param argc  The number of command line arguments of the program
  * \param argv  The contents of the command line arguments of the program
+ * \param printUsage Print a usage string for simulations.
  */
 template <class TypeTag>
 int startLocal(int argc, char **argv,
diff --git a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
index 3c34ff8fd98799c5b4f96cfb1b4b06007946d43a..6a0f7b760e8c8f5baa298a121b559f458e3e52e2 100644
--- a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
+++ b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
@@ -18,7 +18,7 @@
  *****************************************************************************/
 /**
  * \file
- * \brief The problem class for the coupling of an isothermal two-component ZeroEq
+ * \brief The problem which couples an isothermal two-component ZeroEq
  *        and an isothermal two-phase two-component Darcy model.
  */
 #ifndef DUMUX_TWOCZEROEQTWOPTWOCPROBLEM_HH
@@ -95,16 +95,16 @@ SET_TYPE_PROP(TwoCZeroEqTwoPTwoCProblem, LinearSolver, UMFPackBackend<TypeTag>);
 #endif
 }
 
+
 /*!
- * \brief The problem class for the coupling of an isothermal two-component ZeroEq (zeroeq2c)
+ * \ingroup TwoPTwoCZeroEqTwoCModel
+ * \ingroup ImplicitTestProblems
+ * \brief The problem which couples an isothermal two-component ZeroEq (zeroeq2c)
  *        and an isothermal two-phase two-component Darcy model (2p2c).
  *
- *        The problem class for the coupling of a non-isothermal two-component ZeroEq (zeroeq2c)
- *        and an isothermal two-phase two-component Darcy model (2p2c).
- *        It uses the 2p2cCoupling model and the ZeroEq2ccoupling model and provides
- *        the problem specifications for common parameters of the two submodels.
- *        The initial and boundary conditions of the submodels are specified in the two subproblems,
- *        2p2csubproblem.hh and zeroeq2csubproblem.hh, which are accessible via the coupled problem.
+ * It uses the multidomain problem and specifies parameters for the two submodels.
+ * The initial and boundary conditions of the submodels are specified in the two subproblems,
+ * 2p2csubproblem.hh and zeroeq2csubproblem.hh, which are accessible via the coupled problem.
  */
 template <class TypeTag = TTAG(TwoCZeroEqTwoPTwoCProblem) >
 class TwoCZeroEqTwoPTwoCProblem : public MultiDomainProblem<TypeTag>
@@ -169,22 +169,22 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        ElementIterator eendit = mdGrid.template leafend<0>();
-        for (ElementIterator elementIt = mdGrid.template leafbegin<0>();
-             elementIt != eendit; ++elementIt)
+        ElementIterator eEndIt = mdGrid.template leafend<0>();
+        for (ElementIterator eIt = mdGrid.template leafbegin<0>();
+             eIt != eEndIt; ++eIt)
         {
             // this is required for parallelization
             // checks if element is within a partition
-            if (elementIt->partitionType() != Dune::InteriorEntity)
+            if (eIt->partitionType() != Dune::InteriorEntity)
                 continue;
 
-            GlobalPosition globalPos = elementIt->geometry().center();
+            GlobalPosition globalPos = eIt->geometry().center();
 
             if (globalPos[1] > interfacePos_)
-                mdGrid.addToSubDomain(zeroeq2c_,*elementIt);
+                mdGrid.addToSubDomain(zeroeq2c_,*eIt);
             else
                 if(globalPos[0] > noDarcyX_)
-                    mdGrid.addToSubDomain(twoPtwoC_,*elementIt);
+                    mdGrid.addToSubDomain(twoPtwoC_,*eIt);
         }
         mdGrid.preUpdateSubDomains();
         mdGrid.updateSubDomains();
@@ -202,7 +202,7 @@ public:
     void episodeEnd()
     { this->timeManager().startNextEpisode(episodeLength_); }
 
-    //! \copydoc Dumux::CoupledProblem::shouldWriteRestartFile()
+    //! \copydoc Dumux::ImplicitProblem::shouldWriteRestartFile()
     bool shouldWriteRestartFile() const
     {
         return ( ((this->timeManager().timeStepIndex() > 0)
@@ -211,7 +211,7 @@ public:
                 || this->timeManager().episodeWillBeOver());
     }
 
-    //! \copydoc Dumux::CoupledProblem::shouldWriteOutput()
+    //! \copydoc Dumux::ImplicitProblem::shouldWriteOutput()
     bool shouldWriteOutput() const
     {
         return ( ((this->timeManager().timeStepIndex() > 0)
diff --git a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh
index cc3f382933737f0ba45abeeda9ee459e8cab6546..9f8dcb85d1b2853abe6180c22d4556a927297c49 100644
--- a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh
+++ b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cspatialparameters.hh
@@ -58,8 +58,7 @@ SET_TYPE_PROP(TwoCZeroEqTwoPTwoCSpatialParams,
 
 
 /*!
- * \ingroup TwoPTwoCModel
- * \ingroup ZeroEqModel
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \ingroup ImplicitTestProblems
  * \brief Definition of the spatial parameters for
  *        the coupling of an isothermal two-component ZeroEq
diff --git a/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh b/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh
index 7c0ec9394df437ab1a2e56fb7612b21e042f5493..0c34781f6de2e273a0467b8ade45bd40ebbc25e3 100644
--- a/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh
+++ b/test/multidomain/2czeroeq2p2c/2p2csubproblem.hh
@@ -76,14 +76,17 @@ SET_BOOL_PROP(TwoPTwoCSubProblem, ProblemEnableGravity, true);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
+ * \ingroup TwoPTwoCZeroEqTwoCModel
  * \brief Isothermal two-phase two-component porous-medium subproblem
  *        with coupling at the top boundary.
  *
- * \todo update description
+ * The porous-medium subdomain is sized 0.25m times 0.25m. The boundary conditions
+ * are Neumann no-flow everywhere, except at the top, where coupling conditions
+ * are applied to all balance equations. They handle the exchange to the free-flow
+ * subdomain.
  *
- * This sub problem uses the \ref TwoPTwoCModel. It is part of the 2p2c model and
- * is combined with the zeroeq2csubproblem for the free flow domain.
+ * This subproblem uses the \ref TwoPTwoCModel. It is part of a multidomain model and
+ * combined with the zeroeq2csubproblem for the free flow domain.
  */
 template <class TypeTag = TTAG(TwoPTwoCSubProblem) >
 class TwoPTwoCSubProblem : public ImplicitPorousMediaProblem<TypeTag>
@@ -274,11 +277,11 @@ public:
      * \brief Return the initial phase state inside a control volume.
      *
      * \param vertex The vertex
-     * \param globalIdx The index of the global vertex
+     * \param vIdxGlobal The global index of the vertex
      * \param globalPos The global position
      */
     int initialPhasePresence(const Vertex &vertex,
-                             const int &globalIdx,
+                             const int &vIdxGlobal,
                              const GlobalPosition &globalPos) const
     {
         return bothPhases;
diff --git a/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc
index d67b86b6b3eb94789b088de0ad4ecde3b524b297..84e01a85a29f8ff119fc7981c41507f54e693726 100644
--- a/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc
+++ b/test/multidomain/2czeroeq2p2c/test_2czeroeq2p2c.cc
@@ -163,7 +163,7 @@ int startLocal_(int argc, char **argv,
 
     std::string dgfFileName = "temp.dgf";
     Scalar dt, tEnd;
-    Dune::FieldVector<int, dim> nElements;
+    Dune::FieldVector<int, dim> numElements;
     Scalar interfacePos, gradingFactor;
     bool refineTop;
 
@@ -171,8 +171,8 @@ int startLocal_(int argc, char **argv,
         dt = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, DtInitial);
         tEnd = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, TEnd);
 
-        nElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
-        nElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
+        numElements[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsX);
+        numElements[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, CellsY);
         interfacePos = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, InterfacePos);
         gradingFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, GradingY);
         refineTop = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Grid, RefineTop);
@@ -226,7 +226,7 @@ int startLocal_(int argc, char **argv,
     }
 
     Dumux::InterfaceMeshCreator<Grid> interfaceMeshCreator;
-    GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, nElements, interfacePos, gradingFactor, refineTop);
+    GridCreator::gridPtr() = interfaceMeshCreator.create(dgfFileName, numElements, interfacePos, gradingFactor, refineTop);
 
     if (mpiHelper.size() > 1) {
         if (!Dune::Capabilities::isParallel<Grid>::v) {
@@ -299,6 +299,7 @@ int startLocal_(int argc, char **argv,
  *
  * \param argc  The number of command line arguments of the program
  * \param argv  The contents of the command line arguments of the program
+ * \param printUsage  Print a usage string for simulations.
  */
 template <class TypeTag>
 int startLocal(int argc, char **argv,
diff --git a/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh b/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh
index 47471c7a993a79fff91879876fa8bc340d702cea..a7ba88926c295cdd0f5a51cc31bb779cf4e9577a 100644
--- a/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh
+++ b/test/multidomain/2czeroeq2p2c/zeroeq2csubproblem.hh
@@ -88,18 +88,17 @@ SET_SCALAR_PROP(ZeroEq2cSubProblem, FreeFlowSinusTemperaturePeriod, 3600.0);
 
 /*!
  * \ingroup ImplicitTestProblems
- * \ingroup MultidomainProblems
- * \brief ZeroEq2c problem with air flowing from the left to the right.
+ * \ingroup TwoPTwoCZeroEqTwoCModel
+ * \brief Isothermal two-component ZeroEq subproblem with air flowing
+ *        from the left to the right and coupling at the bottom.
  *
- * \todo update test description
- * The zeroeq subdomain is sized 1m times 1m. The boundary conditions for the momentum balances
- * are all set to Dirichlet. The mass balance receives
- * outflow bcs, which are replaced in the localresidual by the sum
- * of the two momentum balances. In the middle of the right boundary,
- * one vertex receives Dirichlet bcs, to set the pressure level.
+ * The free-flow subdomain is sized 0.5m times 0.5m. Dry air is flowing from left (Dirichlet)
+ * to right (outflow), at the right half of the bottom the coupling conditions
+ * are applied to all balance equations. They handle the exchange to the porous-medium
+ * subdomain.
  *
- * This sub problem uses the \ref ZeroEqTwoCModel. It is part of the 2czeroeq2p2c model and
- * is combined with the 2p2csubproblem for the Darcy domain.
+ * This subproblem uses the \ref ZeroEqncModel. It is part of a multidomain model and
+ * combined with the 2p2csubproblem for the porous-medium domain.
  */
 template <class TypeTag>
 class ZeroEq2cSubProblem : public ZeroEqProblem<TypeTag>
@@ -178,11 +177,11 @@ public:
     }
 
     // functions have to be overwritten, otherwise they remain uninitialised
-    //! \copydoc BoxProblem::&bBoxMin()
+    //! \copydoc ImplicitProblem::bBoxMin()
     const GlobalPosition &bBoxMin() const
     { return bBoxMin_; }
 
-    //! \copydoc BoxProblem::&bBoxMax()
+    //! \copydoc ImplicitProblem::bBoxMax()
     const GlobalPosition &bBoxMax() const
     { return bBoxMax_; }