diff --git a/doc/doxygen/modules.txt b/doc/doxygen/modules.txt
index e5585c0d3b4795ab5820af17a560b53df00c6605..aa64dc1d7007913417845e1eb86950f397108e88 100644
--- a/doc/doxygen/modules.txt
+++ b/doc/doxygen/modules.txt
@@ -393,7 +393,7 @@
 /* ***************** Assembly ******************/
 /*!
  * \defgroup Assembly Assembly
- *  TODO: Doc me in modules.txt!
+ * \brief Assembly of linear systems (Jacobian and residual)
  */
 
 /* ***************** Common ******************/
diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index f78804dd4bbcc3418649c6423cdab4c612b1e2b5..5ce3a9f4d012bdabafedef1c67f9851038fe474a 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -18,8 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An assembler for the global linear system for fully implicit models
- *        and cell-centered discretization schemes using Newton's method.
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (box method)
  */
 #ifndef DUMUX_BOX_LOCAL_ASSEMBLER_HH
 #define DUMUX_BOX_LOCAL_ASSEMBLER_HH
@@ -33,16 +34,23 @@
 namespace Dumux {
 
 /*!
- * \ingroup ImplicitModel
- * \brief An assembler for the local contributions (per element) to the global
- *        linear system for fully implicit models and cell-centered discretization schemes.
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (box method)
+ * \tparam TypeTag the TypeTag
+ * \tparam DM the differentiation method to residual compute derivatives
+ * \tparam implicit if to use an implicit or explicit time discretization
  */
 template<class TypeTag,
          DiffMethod DM = DiffMethod::numeric,
          bool implicit = true>
 class BoxLocalAssembler;
 
-
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using numeric differentiation and implicit time discretization
+ */
 template<class TypeTag>
 class BoxLocalAssembler<TypeTag,
                        DiffMethod::numeric,
@@ -442,6 +450,12 @@ private:
 
 }; // implicit BoxAssembler with numeric Jacobian
 
+
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using numeric differentiation and explicit time discretization
+ */
 template<class TypeTag>
 class BoxLocalAssembler<TypeTag,
                        DiffMethod::numeric,
@@ -809,6 +823,12 @@ private:
 
 }; // explicit BoxAssembler with numeric Jacobian
 
+
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using analytic (hard coded) derivatives and implicit time discretization
+ */
 template<class TypeTag>
 class BoxLocalAssembler<TypeTag,
                        DiffMethod::analytic,
@@ -1121,6 +1141,12 @@ private:
 
 }; // implicit BoxAssembler with analytic Jacobian
 
+
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using analytic (hard coded) derivatives and explicit time discretization
+ */
 template<class TypeTag>
 class BoxLocalAssembler<TypeTag,
                        DiffMethod::analytic,
diff --git a/dumux/assembly/boxlocalresidual.hh b/dumux/assembly/boxlocalresidual.hh
index 47db351bb7e7f47f3d7cf586b62c20894932e66d..2dc4423d8f06a0bab5c75e9956ec280141fccecc 100644
--- a/dumux/assembly/boxlocalresidual.hh
+++ b/dumux/assembly/boxlocalresidual.hh
@@ -18,7 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Calculates the residual of models based on the box scheme element-wise.
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Calculates the element-wise residual for the box scheme
  */
 #ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
 #define DUMUX_BOX_LOCAL_RESIDUAL_HH
@@ -29,14 +31,13 @@
 #include <dumux/common/properties.hh>
 #include <dumux/assembly/fvlocalresidual.hh>
 
-namespace Dumux
-{
+namespace Dumux {
+
 /*!
- * \ingroup BoxModel
- * \brief Element-wise calculation of the residual for models
- *        based on the fully implicit box scheme.
- *
- * \todo Please doc me more!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief The element-wise residual for the box scheme
+ * \tparam TypeTag the TypeTag
  */
 template<class TypeTag>
 class BoxLocalResidual : public FVLocalResidual<TypeTag>
@@ -67,6 +68,7 @@ class BoxLocalResidual : public FVLocalResidual<TypeTag>
 public:
     using ParentType::ParentType;
 
+    //! evaluate flux residuals for one sub control volume face and add to residual
     void evalFlux(ElementResidualVector& residual,
                   const Problem& problem,
                   const Element& element,
@@ -91,7 +93,7 @@ public:
         }
     }
 
-
+    //! evaluate flux residuals for one sub control volume face
     ResidualVector evalFlux(const Problem& problem,
                             const Element& element,
                             const FVElementGeometry& fvGeometry,
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 41a6e58c1482d0370da1f97663d5d783f0b907a8..e87077737ecd1e7a1755b6b9598307ecaba896a3 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -18,8 +18,12 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An assembler for the global linear system for fully implicit models
- *        and cell-centered discretization schemes using Newton's method.
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
+ * \tparam TypeTag the TypeTag
+ * \tparam DM the differentiation method to residual compute derivatives
+ * \tparam implicit if to use an implicit or explicit time discretization
  */
 #ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
 #define DUMUX_CC_LOCAL_ASSEMBLER_HH
@@ -34,16 +38,23 @@
 namespace Dumux {
 
 /*!
- * \ingroup ImplicitModel
- * \brief An assembler for the local contributions (per element) to the global
- *        linear system for fully implicit models and cell-centered discretization schemes.
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
+ * \tparam TypeTag the TypeTag
+ * \tparam DM the differentiation method to residual compute derivatives
+ * \tparam implicit if to use an implicit or explicit time discretization
  */
 template<class TypeTag,
          DiffMethod DM = DiffMethod::numeric,
          bool implicit = true>
 class CCLocalAssembler;
 
-
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Cell-centered scheme local assembler using numeric differentiation and implicit time discretization
+ */
 template<class TypeTag>
 class CCLocalAssembler<TypeTag,
                        DiffMethod::numeric,
@@ -586,7 +597,12 @@ private:
     { return gridVolVars.volVars(scv); }
 };
 
-//! Explicit assembler with numeric differentiation
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Cell-centered scheme local assembler using numeric differentiation and explicits time discretization
+ */
 template<class TypeTag>
 class CCLocalAssembler<TypeTag,
                        DiffMethod::numeric,
@@ -926,7 +942,11 @@ private:
     { return gridVolVars.volVars(scv); }
 };
 
-//! implicit assembler using analytic differentiation
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and implicit time discretization
+ */
 template<class TypeTag>
 class CCLocalAssembler<TypeTag,
                        DiffMethod::analytic,
@@ -1187,7 +1207,11 @@ private:
     }
 };
 
-//! explicit assembler using analytic differentiation
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and explicit time discretization
+ */
 template<class TypeTag>
 class CCLocalAssembler<TypeTag,
                        DiffMethod::analytic,
diff --git a/dumux/assembly/cclocalresidual.hh b/dumux/assembly/cclocalresidual.hh
index 15c23241d04cfd7df605a08fd06889ef470e5943..bc13dd2a8961fe643a804fab43c520c7b0f22aa7 100644
--- a/dumux/assembly/cclocalresidual.hh
+++ b/dumux/assembly/cclocalresidual.hh
@@ -18,7 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Calculates the residual of models based on the box scheme element-wise.
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Calculates the element-wise residual for cell-centered discretization schemes
  */
 #ifndef DUMUX_CC_LOCAL_RESIDUAL_HH
 #define DUMUX_CC_LOCAL_RESIDUAL_HH
@@ -28,14 +30,12 @@
 #include <dumux/common/properties.hh>
 #include <dumux/assembly/fvlocalresidual.hh>
 
-namespace Dumux
-{
+namespace Dumux {
+
 /*!
- * \ingroup CCModel
- * \brief Element-wise calculation of the residual for models
- *        based on the fully implicit cell-centered scheme.
- *
- * \todo Please doc me more!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Calculates the element-wise residual for the cell-centered discretization schemes
  */
 template<class TypeTag>
 class CCLocalResidual : public FVLocalResidual<TypeTag>
@@ -54,6 +54,7 @@ class CCLocalResidual : public FVLocalResidual<TypeTag>
 public:
     using ParentType::ParentType;
 
+    //! evaluate the flux residual for a sub control volume face and add to residual
     void evalFlux(ElementResidualVector& residual,
                   const Problem& problem,
                   const Element& element,
@@ -68,6 +69,7 @@ public:
         residual[localScvIdx] += evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
     }
 
+    //! evaluate the flux residual for a sub control volume face
     ResidualVector evalFlux(const Problem& problem,
                             const Element& element,
                             const FVElementGeometry& fvGeometry,
diff --git a/dumux/assembly/diffmethod.hh b/dumux/assembly/diffmethod.hh
index 88e7a17d8ba5356fbfc4fa8930280785f5789628..50b7a8d8bab4afe357226e6d4604eb4556631e2e 100644
--- a/dumux/assembly/diffmethod.hh
+++ b/dumux/assembly/diffmethod.hh
@@ -18,15 +18,21 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An enum class to define various methods available in order to compute
+ * \ingroup Assembly
+ * \brief An enum class to define various differentiation methods available in order to compute
           the derivatives of the residual i.e. the entries in the jacobian matrix.
  */
 #ifndef DUMUX_JACOBIAN_DIFFERENTIATION_METHODS_HH
 #define DUMUX_JACOBIAN_DIFFERENTIATION_METHODS_HH
 
-
 namespace Dumux {
 
+/*!
+ * \ingroup Assembly
+ * \brief Differentiation methods in order to compute the derivatives
+ *        of the residual i.e. the entries in the jacobian matrix.
+ * \todo automatic differentation is not yet implemented
+ */
 enum class DiffMethod
 {
     numeric, analytic, automatic
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index 7d0078cc731f9d5802f3db0b778662349afa3e5a..d8142264a3c4bbc22864c6633c7064cf14340a34 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -18,8 +18,8 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An assembler for the global linear system
- *        for fully implicit models and vertex-centered discretization schemes.
+ * \ingroup Assembly
+ * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
  */
 #ifndef DUMUX_FV_ASSEMBLER_HH
 #define DUMUX_FV_ASSEMBLER_HH
@@ -40,9 +40,11 @@
 namespace Dumux {
 
 /*!
- * \ingroup ImplicitModel
- * \brief An assembler for the global linear system
- *        for fully implicit models and cell-centered discretization schemes.
+ * \ingroup Assembly
+ * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
+ * \tparam TypeTag the TypeTag
+ * \tparam diffMethod the differentiation method to residual compute derivatives
+ * \tparam isImplicit if to use an implicit or explicit time discretization
  */
 template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
 class FVAssembler
diff --git a/dumux/assembly/fvlocalresidual.hh b/dumux/assembly/fvlocalresidual.hh
index bcdc18ebc2fd54c4546485ea94d7d766346d4a44..052b88200aef05891c44eb3a173cdb35a57f3a35 100644
--- a/dumux/assembly/fvlocalresidual.hh
+++ b/dumux/assembly/fvlocalresidual.hh
@@ -18,7 +18,8 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Calculates the element-wise residual of finite-volume models.
+ * \ingroup Assembly
+ * \brief The element-wise residual for finite volume schemes
  */
 #ifndef DUMUX_FV_LOCAL_RESIDUAL_HH
 #define DUMUX_FV_LOCAL_RESIDUAL_HH
@@ -29,15 +30,13 @@
 #include <dumux/common/timeloop.hh>
 #include <dumux/discretization/methods.hh>
 
-namespace Dumux
-{
+namespace Dumux {
 
 /*!
- * \brief Element-wise calculation of the residual matrix for models
- *        using a fully implicit discretization.
- *
+ * \ingroup Assembly
+ * \brief The element-wise residual for finite volume schemes
  * \note This class defines the interface used by the assembler using
- *       static polymorphism.
+ *       static polymorphism. Implementations are specialized for a certain discretization scheme
  */
 template<class TypeTag>
 class FVLocalResidual
@@ -85,8 +84,12 @@ public:
      * This can be used to figure out how much of each conservation
      * quantity is inside the element.
      *
+     * \param problem The problem to solve
      * \param element The DUNE Codim<0> entity for which the storage
      *                term ought to be calculated
+     * \param fvGridGeometry The finite-volume grid geometry
+     * \param gridVariables The grid variables (volume and flux variables)
+     * \param sol The solution vector
      */
     ElementResidualVector evalStorage(const Problem& problem,
                                       const Element &element,
@@ -127,19 +130,17 @@ public:
     /*!
      * \brief Compute the local residual, i.e. the deviation of the
      *        equations from zero for instationary problems.
+     *
      * \param problem The problem to solve
-
      * \param element The DUNE Codim<0> entity for which the residual
      *                ought to be calculated
      * \param fvGeometry The finite-volume geometry of the element
-     * \param prevVolVars The volume averaged variables for all
-     *                    sub-control volumes of the element at the previous
-     *                    time level
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
-     * \param bcTypes The types of the boundary conditions for all
-     *                vertices of the element
+     * \param prevElemVolVars The volume averaged variables for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param bcTypes The types of the boundary conditions for all boundary entities of an element
+     * \param elemFluxVarsCache The flux variable caches for the element stencil
      */
     ElementResidualVector eval(const Problem& problem,
                                const Element& element,
@@ -176,19 +177,17 @@ public:
     /*!
      * \brief Compute the storage local residual, i.e. the deviation of the
      *        storage term from zero for instationary problems.
+     *
      * \param problem The problem to solve
-
      * \param element The DUNE Codim<0> entity for which the residual
      *                ought to be calculated
      * \param fvGeometry The finite-volume geometry of the element
-     * \param prevVolVars The volume averaged variables for all
-     *                    sub-control volumes of the element at the previous
-     *                    time level
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
-     * \param bcTypes The types of the boundary conditions for all
-     *                vertices of the element
+     * \param prevElemVolVars The volume averaged variables for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param bcTypes The types of the boundary conditions for all boundary entities of an element
+     * \param elemFluxVarsCache The flux variable caches for the element stencil
      */
     ElementResidualVector evalStorage(const Problem& problem,
                                       const Element& element,
@@ -218,16 +217,15 @@ public:
     /*!
      * \brief Compute the local residual, i.e. the deviation of the
      *        equations from zero for stationary problem.
+     *
      * \param problem The problem to solve
-
      * \param element The DUNE Codim<0> entity for which the residual
      *                ought to be calculated
      * \param fvGeometry The finite-volume geometry of the element
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
-     * \param bcTypes The types of the boundary conditions for all
-     *                vertices of the element
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param bcTypes The types of the boundary conditions for all boundary entities of an element
+     * \param elemFluxVarsCache The flux variable caches for the element stencil
      */
     ElementResidualVector eval(const Problem& problem,
                                const Element& element,
@@ -261,14 +259,16 @@ public:
 
     /*!
      * \name Model specific interface
-     * \note The following method are the model specific implementations of the actual residual
+     * \note The following method are the model specific implementations of the local residual
      */
     // \{
 
     /*!
      * \brief Calculate the source term of the equation
      *
-     * \param scv The sub-control volume over which we integrate the source term
+     * \param problem The problem to solve
+     * \param scv The sub-control volume over which we integrate the storage term
+     * \param volVars The volume variables associated with the scv
      * \note has to be implemented by the model specific residual class
      *
      */
@@ -282,6 +282,11 @@ public:
     /*!
      * \brief Calculate the source term of the equation
      *
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param elemVolVars The volume variables associated with the element stencil
      * \param scv The sub-control volume over which we integrate the source term
      * \note This is the default implementation for all models as sources are computed
      *       in the user interface of the problem
@@ -307,7 +312,14 @@ public:
     /*!
      * \brief Calculate the source term of the equation
      *
-     * \param scv The sub-control volume over which we integrate the source term
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param elemVolVars The volume variables associated with the element stencil
+     * \param scvf The sub-control volume over which we integrate the flux
+     * \param elemFluxVarsCache the flux variable caches for the element's flux stencils
+     *
      * \note has to be implemented by the model specific residual class
      *
      */
@@ -329,6 +341,21 @@ public:
      */
     // \{
 
+    /*!
+     * \brief Compute the storage local residual, i.e. the deviation of the
+     *        storage term from zero for instationary problems.
+     *
+     * \param residual The residual vector to fill
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param prevElemVolVars The volume averaged variables for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param scv The sub control volume the storage term is integrated over
+     */
     void evalStorage(ElementResidualVector& residual,
                      const Problem& problem,
                      const Element& element,
@@ -344,7 +371,7 @@ public:
         // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
         // euler as time discretization.
         //
-        // We might need a more explicit way for
+        // TODO: We might need a more explicit way for
         // doing the time discretization...
 
         //! Compute storage with the model specific storage residual
@@ -361,6 +388,19 @@ public:
         residual[scv.indexInElement()] += storage;
     }
 
+    /*!
+     * \brief Compute the source local residual, i.e. the deviation of the
+     *        source term from zero.
+     *
+     * \param residual The residual vector to fill
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param scv The sub control volume the source term is integrated over
+     */
     void evalSource(ElementResidualVector& residual,
                     const Problem& problem,
                     const Element& element,
@@ -377,6 +417,22 @@ public:
         residual[scv.indexInElement()] -= source;
     }
 
+    /*!
+     * \brief Compute the source local residual, i.e. the deviation of the
+     *        source term from zero and add to the residual
+     * \note This is implemented for different discretization schemes
+     *
+     * \param residual The residual vector to fill
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param elemVolVars The volume averaged variables for all
+     *                    sub-control volumes of the element at the current  time level
+     * \param elemBcTypes the boundary types for the boundary entities of an elements
+     * \param elemFluxVarsCache The flux variable caches for the element stencil
+     * \param scvf The sub control volume face the flux term is integrated over
+     */
     void evalFlux(ElementResidualVector& residual,
                   const Problem& problem,
                   const Element& element,
@@ -386,6 +442,19 @@ public:
                   const ElementFluxVariablesCache& elemFluxVarsCache,
                   const SubControlVolumeFace& scvf) const {}
 
+    /*!
+     * \brief Compute the flux local residual, i.e. the deviation of the
+     *        flux term from zero.
+     *
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param elemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param elemFluxVarsCache The flux variable caches for the element stencil
+     * \param scvf The sub control volume face the flux term is integrated over
+     */
     ResidualVector evalFlux(const Problem& problem,
                             const Element& element,
                             const FVElementGeometry& fvGeometry,
@@ -396,6 +465,14 @@ public:
         return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
     }
 
+    //\}
+
+    /*!
+     * \name Interfaces for analytic Jacobian computation
+     */
+    // \{
+
+    //! Compute the derivative of the storage residual
     template<class PartialDerivativeMatrix>
     void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
                                const Problem& problem,
@@ -407,6 +484,7 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic storage derivative");
     }
 
+    //! Compute the derivative of the source residual
     template<class PartialDerivativeMatrix>
     void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
                               const Problem& problem,
@@ -418,6 +496,7 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic source derivative");
     }
 
+    //! Compute the derivative of the flux residual
     template<class PartialDerivativeMatrices, class T = TypeTag>
     std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) != DiscretizationMethods::Box, void>
     addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
@@ -431,6 +510,7 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for cell-centered models");
     }
 
+    //! Compute the derivative of the flux residual for the box method
     template<class JacobianMatrix, class T = TypeTag>
     std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Box, void>
     addFluxDerivatives(JacobianMatrix& A,
@@ -444,6 +524,7 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for box models");
     }
 
+    //! Compute the derivative of the Dirichlet flux residual for cell-centered schemes
     template<class PartialDerivativeMatrices>
     void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                      const Problem& problem,
@@ -456,6 +537,7 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic Dirichlet flux derivative");
     }
 
+    //! Compute the derivative of Robin type boundary conditions ("solution dependent Neumann")
     template<class PartialDerivativeMatrices>
     void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                  const Problem& problem,
@@ -468,15 +550,17 @@ public:
         DUNE_THROW(Dune::NotImplemented, "analytic Robin flux derivative");
     }
 
+    //\}
+
     /*!
-     * \brief Sets the solution from which to start the time integration. Has to be
-     *        called prior to assembly for time-dependent problems.
+     * \brief Sets the solution from which to start the time integration.
+     * \note Has to be called prior to assembly for time-dependent problems.
      */
     void setPreviousSolution(const SolutionVector& u)
     { prevSol_ = &u; }
 
     /*!
-     * \brief Return the solution that has been set as the previous one.
+     * \brief Return the solution of the previous time level that has been set
      */
     const SolutionVector& prevSol() const
     {
@@ -485,16 +569,21 @@ public:
     }
 
     /*!
-     * \brief If no solution has been set, we treat the problem as stationary.
+     * \brief if the residual is stationary (storage residual is strictly zero)
+     * \note If no solution has been set, we treat the problem as stationary.
      */
     bool isStationary() const
     { return !prevSol_; }
 
     // \}
 protected:
+    //! the timeloop for instationary problems
+    //! calling this for stationary leads to undefined behaviour
     TimeLoop& timeLoop()
     { return *timeLoop_; }
 
+    //! the timeloop for instationary problems
+    //! calling this for stationary leads to undefined behaviour
     const TimeLoop& timeLoop() const
     { return *timeLoop_; }
 
@@ -505,8 +594,8 @@ protected:
     { return *static_cast<const Implementation*>(this); }
 
 private:
-    std::shared_ptr<TimeLoop> timeLoop_;
-    const SolutionVector* prevSol_;
+    std::shared_ptr<TimeLoop> timeLoop_; //!< the timeloop for instationary problems
+    const SolutionVector* prevSol_; //!< the solution of the previous time level for instationary problems
 };
 
 } // end namespace Dumux
diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
index 7afc878ba075cf2199f323d8b3dcbb75393b7acb..7042b115c2a7e310ff76590de67a02116d3cc2b3 100644
--- a/dumux/assembly/staggeredfvassembler.hh
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -18,8 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An assembler for the global linear system
- *        for fully implicit models and vertex-centered discretization schemes.
+ * \ingroup Assembly
+ * \ingroup StaggeredDiscretization
+ * \brief A linear system assembler (residual and Jacobian) for staggered finite volume schemes
  */
 #ifndef DUMUX_STAGGERED_FV_ASSEMBLER_HH
 #define DUMUX_STAGGERED_FV_ASSEMBLER_HH
@@ -39,9 +40,12 @@
 namespace Dumux {
 
 /*!
- * \ingroup ImplicitModel
- * \brief An assembler for the global linear system
- *        for fully implicit models and cell-centered discretization schemes.
+ * \ingroup Assembly
+ * \ingroup StaggeredDiscretization
+ * \brief A linear system assembler (residual and Jacobian) for staggered finite volume schemes
+ * \tparam TypeTag the TypeTag
+ * \tparam diffMethod the differentiation method to residual compute derivatives
+ * \tparam isImplicit if to use an implicit or explicit time discretization
  */
 template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
 class StaggeredFVAssembler
@@ -255,7 +259,6 @@ public:
     /*!
      * \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
      */
-
     void setJacobianPattern()
     {
         // resize the jacobian and the residual
@@ -360,14 +363,14 @@ public:
     { return localResidual_; }
 
 private:
-    // reset the residual to 0.0
+    //! reset the residual to 0.0
     void resetResidual_()
     {
         (*residual_)[cellCenterIdx] = 0.0;
         (*residual_)[faceIdx] = 0.0;
     }
 
-    // reset the jacobian to 0.0
+    //! reset the jacobian to 0.0
     void resetJacobian_()
     {
         (*jacobian_)[cellCenterIdx][cellCenterIdx] = 0.0;
@@ -376,23 +379,23 @@ private:
         (*jacobian_)[faceIdx][faceIdx] = 0.0;
     }
 
-    // pointer to the problem to be solved
+    //! pointer to the problem to be solved
     std::shared_ptr<const Problem> problem_;
 
-    // the finite volume geometry of the grid
+    //! the finite volume geometry of the grid
     std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
 
-    // the variables container for the grid
+    //! the variables container for the grid
     std::shared_ptr<GridVariables> gridVariables_;
 
-    // shared pointers to the jacobian matrix and residual
+    //! shared pointers to the jacobian matrix and residual
     std::shared_ptr<JacobianMatrix> jacobian_;
     std::shared_ptr<SolutionVector> residual_;
 
-    // class computing the residual of an element
+    //! class computing the residual of an element
     LocalResidual localResidual_;
 
-    // if this assembler is assembling a time dependent problem
+    //! if this assembler is assembling a time dependent problem
     bool stationary_;
 };
 
diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
index fbd52f3d4222d2433ae1d48a56025f9df292659f..1a95b101cd1ddc6774ce39e34151bce8f724d79b 100644
--- a/dumux/assembly/staggeredlocalassembler.hh
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -18,8 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief An assembler for the global linear system for fully implicit models
- *        and cell-centered discretization schemes using Newton's method.
+ * \ingroup StaggeredDiscretization
+ * \ingroup Assembly
+ * \brief An assembler for Jacobian and residual contribution per element (staggered FV method)
  */
 #ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
 #define DUMUX_CC_LOCAL_ASSEMBLER_HH
@@ -35,24 +36,30 @@
 #include <dune/common/version.hh>
 
 #if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
-    #include <dune/common/hybridutilities.hh>
+#include <dune/common/hybridutilities.hh>
 #else
-    #include <dumux/common/intrange.hh>
+#include <dumux/common/intrange.hh>
 #endif
 
 namespace Dumux {
 
 /*!
- * \ingroup ImplicitModel
- * \brief An assembler for the local contributions (per element) to the global
- *        linear system for fully implicit models and cell-centered discretization schemes.
+ * \ingroup Assembly
+ * \ingroup StaggeredDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (staggered FV method)
+ * \tparam TypeTag the TypeTag
+ * \tparam DM the differentiation method to residual compute derivatives
+ * \tparam implicit if to use an implicit or explicit time discretization
  */
 template<class TypeTag,
          DiffMethod DM = DiffMethod::numeric,
          bool implicit = true>
 class StaggeredLocalAssembler;
 
-
+/*!
+ * \ingroup Assembly
+ * \brief Staggerd FV scheme local assembler using numeric differentiation and implicit time discretization
+ */
 template<class TypeTag>
 class StaggeredLocalAssembler<TypeTag,
                        DiffMethod::numeric,
@@ -548,11 +555,14 @@ protected:
 
     /*!
      * \brief Computes the epsilon used for numeric differentiation
-     *        for a given value of a primary variable.
+     *        for a given value of a primary variable for
+     *        an equation residual with respect to a primary variable
      *
      * \param priVar The value of the primary variable
+     * \param eqIdx The equation index of the equation
+     * \param pvIdx The index of the primary variable
      */
-    static Scalar numericEpsilon(const Scalar priVar, const int idx1, const int idx2)
+    static Scalar numericEpsilon(const Scalar priVar, const int eqIdx, const int pvIdx)
     {
         // define the base epsilon as the geometric mean of 1 and the
         // resolution of the scalar type. E.g. for standard 64 bit
@@ -566,7 +576,7 @@ protected:
         const std::array<std::array<Scalar, 2>, 2> baseEps_ = BaseEpsilon::getEps();
 
 
-        static const Scalar baseEps = baseEps_[idx1][idx2];
+        static const Scalar baseEps = baseEps_[eqIdx][pvIdx];
         assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
         // the epsilon value used for the numeric differentiation is
         // now scaled by the absolute value of the primary variable...
diff --git a/dumux/assembly/staggeredlocalresidual.hh b/dumux/assembly/staggeredlocalresidual.hh
index 5c0179a669ffa64cefde07b5a3984395777a34df..2f14bf6de631df5690241eadb7cac2fd643bfc05 100644
--- a/dumux/assembly/staggeredlocalresidual.hh
+++ b/dumux/assembly/staggeredlocalresidual.hh
@@ -18,7 +18,9 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Calculates the residual of models based on the box scheme element-wise.
+ * \ingroup StaggeredDiscretization
+ * \ingroup Assembly
+ * \brief Calculates the element-wise residual for the staggered FV scheme
  */
 #ifndef DUMUX_STAGGERED_LOCAL_RESIDUAL_HH
 #define DUMUX_STAGGERED_LOCAL_RESIDUAL_HH
@@ -26,20 +28,16 @@
 #include <dumux/common/valgrind.hh>
 #include <dumux/common/timeloop.hh>
 
-namespace Dumux
-{
+namespace Dumux {
 
 namespace Properties
 {
     NEW_PROP_TAG(ElementFaceVariables);
 }
 /*!
- * \ingroup CCModel
- * \ingroup StaggeredLocalResidual
- * \brief Element-wise calculation of the residual for models
- *        based on the fully implicit cell-centered scheme.
- *
- * \todo Please doc me more!
+ * \ingroup StaggeredDiscretization
+ * \ingroup Assembly
+ * \brief Calculates the element-wise residual for the staggered FV scheme
  */
 template<class TypeTag>
 class StaggeredLocalResidual
@@ -93,28 +91,32 @@ public:
     , prevSol_(nullptr)
     {}
 
-     /*!
+    /*!
      * \name User interface
      * \note The following methods are usually expensive to evaluate
      *       They are useful for outputting residual information.
      */
     // \{
 
-     /*!
-     * \brief Compute the local residual, i.e. the deviation of the
+    /*!
+     * \brief Compute the local cell residual, i.e. the deviation of the
      *        equations from zero for a transient problem.
      *
+     * \param problem The problem to solve
      * \param element The DUNE Codim<0> entity for which the residual
      *                ought to be calculated
      * \param fvGeometry The finite-volume geometry of the element
-     * \param prevVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the previous
-     *                   time level
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
+     * \param prevElemVolVars The volume averaged variables (in the cell) for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables (in the cell) for all
+     *                       sub-control volumes of the element at the current time level
+     * \param prevElemFaceVars The volume averaged variables (on the face) for all
+     *                         sub-control volumes of the element at the previous time level
+     * \param curElemFaceVars The volume averaged variables (on the face) for all
+     *                        sub-control volumes of the element at the current time level
      * \param bcTypes The types of the boundary conditions for all
      *                vertices of the element
+     * \param elemFluxVarsCache The cache of the flux variables
      */
     auto evalCellCenter(const Problem& problem,
                         const Element &element,
@@ -142,21 +144,26 @@ public:
         return residual;
     }
 
-     /*!
-     * \brief Compute the local residual, i.e. the deviation of the
+    /*!
+     * \brief Compute the local face residual, i.e. the deviation of the
      *        equations from zero for a transient problem.
      *
+     * \param problem The problem to solve
      * \param element The DUNE Codim<0> entity for which the residual
      *                ought to be calculated
      * \param fvGeometry The finite-volume geometry of the element
-     * \param prevVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the previous
-     *                   time level
-     * \param curVolVars The volume averaged variables for all
-     *                   sub-control volumes of the element at the current
-     *                   time level
+     * \param scvf The sub control volume face
+     * \param prevElemVolVars The volume averaged variables (in the cell) for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables (in the cell) for all
+     *                       sub-control volumes of the element at the current time level
+     * \param prevElemFaceVars The volume averaged variables (on the face) for all
+     *                         sub-control volumes of the element at the previous time level
+     * \param curElemFaceVars The volume averaged variables (on the face) for all
+     *                        sub-control volumes of the element at the current time level
      * \param bcTypes The types of the boundary conditions for all
      *                vertices of the element
+     * \param elemFluxVarsCache The cache of the flux variables
      */
     auto evalFace(const Problem& problem,
                   const Element &element,