diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index a2427e9813d569f9a75397ca94b59233736431e5..864299844e0da459d83c400eb42f78ac028f6468 100644
--- a/dumux/freeflow/stokes/stokesfluxvariables.hh
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -43,11 +43,11 @@ namespace Dumux
  * \ingroup BoxStokesModel
  * \ingroup BoxFluxVariables
  * \brief This template class contains the data which is required to
- *        calculate the fluxes of the fluid phase over a face of a
- *        finite volume for the Stokes box model.
+ *        calculate the mass and momentum fluxes over the face of a
+ *        sub-control volume for the Stokes box model.
  *
- * This means pressure gradients, phase densities and viscosities, etc. at
- * the integration point of the sub-control-volume face
+ * This means pressure gradients, phase densities, viscosities, etc.
+ * at the integration point of the sub-control-volume face
  */
 template <class TypeTag>
 class StokesFluxVariables
@@ -173,7 +173,7 @@ public:
 
     /*!
      * \brief Return the average volume of the upstream and the downstream sub-control volume;
-     *        this is required for the stabilization
+     *        this is required for the stabilization.
      */
     const Scalar averageSCVVolume() const
     {
@@ -182,21 +182,21 @@ public:
     }
 
     /*!
-     * \brief Return pressure \f$\mathrm{[Pa]}\f$ at the integration
+     * \brief Return the pressure \f$\mathrm{[Pa]}\f$ at the integration
      *        point.
      */
     Scalar pressureAtIP() const
     { return pressureAtIP_; }
 
     /*!
-     * \brief Return density \f$\mathrm{[kg/m^3]}\f$ at the integration
+     * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ at the integration
      *        point.
      */
     Scalar densityAtIP() const
     { return densityAtIP_; }
 
     /*!
-     * \brief Return viscosity \f$\mathrm{[m^2/s]}\f$ at the integration
+     * \brief Return the viscosity \f$\mathrm{[m^2/s]}\f$ at the integration
      *        point.
      */
     Scalar viscosityAtIP() const
@@ -210,47 +210,47 @@ public:
     { return normalVelocityAtIP_; }
 
     /*!
-     * \brief Return the pressure gradient at the integration
-     *        point.
+     * \brief Return the pressure gradient at the integration point.
      */
     const ScalarGradient &pressureGradAtIP() const
     { return pressureGradAtIP_; }
 
     /*!
-     * \brief Return the velocity at the integration
-     *        point.
+     * \brief Return the velocity vector at the integration point.
      */
     const VelocityVector &velocityAtIP() const
     { return velocityAtIP_; }
 
     /*!
      * \brief Return the velocity gradient at the integration
-     *        point.
+     *        point of a face.
      */
     const VectorGradient &velocityGradAtIP() const
     { return velocityGradAtIP_; }
 
     /*!
-     * \brief Return the divergence of the normal velocity at the integration
-     *        point.
+     * \brief Return the divergence of the normal velocity at the
+     *        integration point.
      */
     Scalar velocityDivAtIP() const
     { return velocityDivAtIP_; }
 
     /*!
-     * \brief Return the local index of the upstream control volume
-     *        for a given phase.
+     * \brief Return the local index of the upstream sub-control volume.
      */
     int upstreamIdx() const
     { return upstreamIdx_; }
 
     /*!
-     * \brief Return the local index of the downstream control volume
-     *        for a given phase.
+     * \brief Return the local index of the downstream sub-control volume.
      */
     int downstreamIdx() const
     { return downstreamIdx_; }
 
+    /*!
+     * \brief Indicates if a face is on a boundary. Used for in the
+     *        face() method (e.g. for outflow boundary conditions).
+     */
     bool onBoundary() const
     { return onBoundary_; }
 
diff --git a/dumux/freeflow/stokes/stokesindices.hh b/dumux/freeflow/stokes/stokesindices.hh
index 7d0c74bc0af632d26053c00817785d557e9eae58..3afa6bd193dd86698df84b8518ee0d9c7508b8fc 100644
--- a/dumux/freeflow/stokes/stokesindices.hh
+++ b/dumux/freeflow/stokes/stokesindices.hh
@@ -24,7 +24,7 @@
 /*!
  * \file
  *
- * \brief Defines the indices required for the Stokes BOX model.
+ * \brief Defines the indices required for the Stokes box model.
  */
 #ifndef DUMUX_STOKES_INDICES_HH
 #define DUMUX_STOKES_INDICES_HH
diff --git a/dumux/freeflow/stokes/stokeslocaljacobian.hh b/dumux/freeflow/stokes/stokeslocaljacobian.hh
index 841d6a08c7339fa21bf311af37974950e989b2ea..f934c83756f472a02da9e45dcd5c261daf48dba8 100644
--- a/dumux/freeflow/stokes/stokeslocaljacobian.hh
+++ b/dumux/freeflow/stokes/stokeslocaljacobian.hh
@@ -38,12 +38,14 @@ namespace Dumux
  * \brief Element-wise calculation of the jacobian matrix for models
  *        based on the box scheme.
  *
- * This overloads the numericEpsilon method of the boxlocaljacobian
+ * This overloads the numericEpsilon method of the boxlocaljacobian.
+ * The momentum balance equation uses larger epsilons than the rest.
  */
 template<class TypeTag>
 class StokesLocalJacobian : public BoxLocalJacobian<TypeTag>
 {
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
 
 public:
     //! \copydoc BoxLocalJacobian::numericEpsilon()
@@ -51,7 +53,7 @@ public:
                           int pvIdx) const
     {
         Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
-        if (pvIdx == 0 || pvIdx == 1){
+        if (pvIdx < GridView::dimension){
             return 1e-7*(std::abs(pv) + 1);
         }
         return 1e-9*(std::abs(pv) + 1);
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index 7056b1a34f090ebeac85edf23f8ac4d7ba8bdee1..5a983e145ec9081d398a96a5e6f4a4166abb8a99 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -24,7 +24,7 @@
  * \file
  *
  * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the stokes box model.
+ *        using the Stokes box model.
  */
 
 #ifndef DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
@@ -52,7 +52,7 @@ namespace Dumux
  *        using the Stokes box model.
  *
  * This class is also used for the non-isothermal and the two-component Stokes
- * model, which means that it uses static polymorphism.
+ * model (static polymorphism).
  */
 template<class TypeTag>
 class StokesLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
@@ -84,7 +84,6 @@ protected:
     typedef Dune::FieldVector<Scalar, dim> FieldVector;
 
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
@@ -108,7 +107,7 @@ protected:
     };
 
     /*!
-     * \brief Evaluate the amount all conservation quantities
+     * \brief Evaluates the amount of all conservation quantities
      *        (mass and momentum) within a sub-control volume.
      *
      * The result should be averaged over the volume (e.g. phase mass
@@ -149,7 +148,7 @@ protected:
      * \param flux The flux over the SCV (sub-control-volume) face
      * \param faceIdx The index of the SCV face (may also be a boundary face)
      * \param onBoundary Indicates, if the flux is evaluated on a boundary face. If it is true,
-     *        the created fluxVars object cotains boundary variables evaluated at the IP of the
+     *        the created fluxVars object contains boundary variables evaluated at the IP of the
      *        boundary face
      */
     void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const
@@ -173,7 +172,7 @@ protected:
      *        a face of a sub-control volume.
      *
      * \param flux The advective flux over the sub-control-volume face for each component
-     * \param fluxVars The flux variables at the current SCV face
+     * \param fluxVars The flux variables at the current SCV/boundary face
      */
     void computeAdvectiveFlux(PrimaryVariables &flux,
                               const FluxVariables &fluxVars) const
@@ -250,14 +249,14 @@ protected:
 
     /*!
      * \brief Adds the diffusive flux to the flux vector over
-     *        the face of a sub-control volume.
+     *        a SCV face or a boundary face.
      *
      * It doesn't do anything in the Stokes model but is used by the
      * transport and non-isothermal models to calculate diffusive and
-     * conductive fluxes
+     * conductive fluxes.
      *
-     * \param flux The diffusive flux over the sub-control-volume face for each component
-     * \param fluxVars The flux variables at the current sub control volume face
+     * \param flux The diffusive flux over the SCV face or boundary face for each component
+     * \param fluxVars The flux variables at the current SCV/boundary face
      */
     void computeDiffusiveFlux(PrimaryVariables &flux,
                               const FluxVariables &fluxVars) const
@@ -269,7 +268,7 @@ protected:
      *        and the gravity term evaluated.
      *
      * \param q The source/sink in the sub control volume for each component
-     * \param localVertexIdx The index of the sub-control volume
+     * \param localVertexIdx The local index of the sub-control volume
      */
     void computeSource(PrimaryVariables &q, int localVertexIdx)
     {
@@ -434,7 +433,7 @@ protected:
             // call evalNeumannSegment_() of the base class first
             ParentType::evalNeumannSegment_(isIt, scvIdx, boundaryFaceIdx);
 
-            // temporary vector to store the neumann boundary fluxes
+            // temporary vector to store the Neumann boundary fluxes
             PrimaryVariables values(0.0);
             if (momentumBalanceHasNeumann_(bcTypes))
             {
@@ -477,8 +476,7 @@ protected:
     }
 
     /*!
-     * \brief Evaluate outflow boundary conditions for a single sub-control
-     *        volume face.
+     * \brief Evaluate outflow boundary conditions for a single SCV face on the boundary.
      */
     void evalOutflowSegment_(const IntersectionIterator &isIt,
                              const int scvIdx,
@@ -576,7 +574,7 @@ protected:
 
     /*!
      * \brief Interpolate the pressure at corner points of the grid, thus taking the degree of freedom there. 
-     * 		  This is required due to stability reasons.s
+     * 		  This is required due to stability reasons.
      */
     void interpolateCornerPoints_(const BoundaryTypes &bcTypes, const int vertexIdx)
     {
@@ -601,7 +599,7 @@ protected:
 
     /*!
      * \brief Replace the local residual of the mass balance equation by
-     *        the sum of the residuals of the momentum balance equation
+     *        the sum of the residuals of the momentum balance equation.
      */
     void replaceMassbalanceResidual_(const FieldVector& momentumResidual,
                                      FieldVector& averagedNormal,
@@ -617,7 +615,7 @@ protected:
 
     /*!
      * \brief Returns true, if all boundary conditions for the momentum balance
-     *        at the considered vertex are dirichlet
+     *        at the considered vertex are Dirichlet.
      */
     bool momentumBalanceDirichlet_(const BoundaryTypes& bcTypes) const
     {
@@ -628,7 +626,7 @@ protected:
     }
 
     /*!
-     * \brief Returns true, if at least one boundary condition of the momentum balance is neumann
+     * \brief Returns true, if at least one boundary condition of the momentum balance is Neumann.
      */
     bool momentumBalanceHasNeumann_(const BoundaryTypes& bcTypes) const
     {
@@ -639,7 +637,7 @@ protected:
     }
 
     /*!
-     * \brief Returns true, if all boundary conditions for the momentum balance are outlow
+     * \brief Returns true, if all boundary conditions for the momentum balance are outflow.
      */
     bool momentumBalanceOutflow_(const BoundaryTypes& bcTypes) const
     {
diff --git a/dumux/freeflow/stokes/stokesmodel.hh b/dumux/freeflow/stokes/stokesmodel.hh
index c5e0269cf15595927b3beee46a9b14508b911748..a4f64abdc3dc126b12a4f96e09b0740ebe588e25 100644
--- a/dumux/freeflow/stokes/stokesmodel.hh
+++ b/dumux/freeflow/stokes/stokesmodel.hh
@@ -28,6 +28,8 @@
 #ifndef DUMUX_STOKES_MODEL_HH
 #define DUMUX_STOKES_MODEL_HH
 
+#include <dumux/boxmodels/common/boxmodel.hh>
+
 #include "stokeslocalresidual.hh"
 #include "stokesnewtoncontroller.hh"
 #include "stokeslocaljacobian.hh"
@@ -38,7 +40,7 @@ namespace Dumux
 {
 /*!
  * \ingroup BoxStokesModel
- * \brief Adaption of the box scheme to the stokes model.
+ * \brief Adaption of the box scheme to the Stokes model.
  *
  * This model implements laminar Stokes flow of a single fluid, solving a momentum balance:
  * \f[
@@ -67,9 +69,9 @@ class StokesModel : public BoxModel<TypeTag>
         dim = GridView::dimension,
         dimWorld = GridView::dimensionworld
     };
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
 
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
 
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
@@ -84,10 +86,14 @@ public:
     /*!
      * \brief Calculate the fluxes across a certain layer in the domain.
      * The layer is situated perpendicular to the coordinate axis "coord" and cuts
-     * the axis at the value "coordValue"
+     * the axis at the value "coordVal".
      *
+     * \param globalSol The global solution vector
+     * \param flux A vector to store the flux
+     * \param axis The dimension, perpendicular to which the layer is situated
+     * \param coordVal The (Scalar) coordinate on the axis, at which the layer is situated
      */
-    void calculateFluxAcrossLayer (const SolutionVector &globalSol, Dune::FieldVector<Scalar, 2> &flux, int coord, Scalar coordVal)
+    void calculateFluxAcrossLayer(const SolutionVector &globalSol, Dune::FieldVector<Scalar, numEq> &flux, int axis, Scalar coordVal)
     {
         GlobalPosition globalI, globalJ;
         PrimaryVariables tmpFlux(0.0);
@@ -112,9 +118,9 @@ public:
             bool hasRight = false;
             for (int i = 0; i < fvElemGeom.numVertices; i++) {
                 const GlobalPosition &global = fvElemGeom.subContVol[i].global;
-                if (globalI[coord] < coordVal)
+                if (globalI[axis] < coordVal)
                     hasLeft = true;
-                else if (globalI[coord] >= coordVal)
+                else if (globalI[axis] >= coordVal)
                     hasRight = true;
             }
             if (!hasLeft || !hasRight)
@@ -124,7 +130,7 @@ public:
                 const int &globalIdx =
                     this->vertexMapper().map(*elemIt, i, dim);
                 const GlobalPosition &global = fvElemGeom.subContVol[i].global;
-                if (globalI[coord] < coordVal)
+                if (globalI[axis] < coordVal)
                     flux += this->localResidual().residual(i);
             }
         }
@@ -132,11 +138,7 @@ public:
         flux = this->problem_.gridView().comm().sum(flux);
     }
 
-    /*!
-     * \brief Calculate mass in the whole model domain
-     *         and get minimum and maximum values of primary variables
-     *
-     */
+    DUMUX_DEPRECATED_MSG("outdated function; use globalStorage() of the model")
     void calculateMass(const SolutionVector &sol, Dune::FieldVector<Scalar, 2> &mass)
     {
         //         const DofMapper &dofMapper = this->dofEntityMapper();
diff --git a/dumux/freeflow/stokes/stokesnewtoncontroller.hh b/dumux/freeflow/stokes/stokesnewtoncontroller.hh
index 4b7e46a31d390414d747646408e629ad62b4255c..fd78e2d336e876911f917d5dc93b84aed4fcaa80 100644
--- a/dumux/freeflow/stokes/stokesnewtoncontroller.hh
+++ b/dumux/freeflow/stokes/stokesnewtoncontroller.hh
@@ -33,7 +33,7 @@ namespace Dumux {
 /*!
  * \ingroup BoxStokesModel
  * \ingroup Newton
- * \brief A Stokes-specific controller for the newton solver, which sets
+ * \brief A Stokes-specific controller for the nonlinear Newton solver, which sets
  *        different parameters for the relative tolerance, target steps and maximum steps.
  */
 template <class TypeTag>
diff --git a/dumux/freeflow/stokes/stokesproblem.hh b/dumux/freeflow/stokes/stokesproblem.hh
index 8dbda9354f0f3450f38322637ded151efd876dc8..0f86e23a1c889310db4301f05f033bdd4c2afe8b 100644
--- a/dumux/freeflow/stokes/stokesproblem.hh
+++ b/dumux/freeflow/stokes/stokesproblem.hh
@@ -22,7 +22,7 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Base class for all stokes problems which use the box scheme
+ * \brief Base class for all stokes problems which use the box scheme.
  */
 #ifndef DUMUX_STOKES_PROBLEM_HH
 #define DUMUX_STOKES_PROBLEM_HH
@@ -35,9 +35,9 @@ namespace Dumux
 {
 /*!
  * \ingroup BoxStokesProblems
- * \brief Base class for all problems which use the stokes box model
+ * \brief Base class for all problems which use the Stokes box model.
  *
- * \todo Please doc me more!
+ * This implements gravity (if desired) and a function returning the temperature.
  */
 template<class TypeTag>
 class StokesProblem : public BoxProblem<TypeTag>
@@ -97,7 +97,7 @@ public:
      * \brief Evaluate the intrinsic permeability
      *        at the corner of a given element
      *
-     * \return permeability in x-direction
+     * \return (Scalar) permeability
      */
     Scalar permeability(const Element &element,
                         const FVElementGeometry &fvElemGeom,
diff --git a/dumux/freeflow/stokes/stokesproperties.hh b/dumux/freeflow/stokes/stokesproperties.hh
index daab9762d064073c07cbbbb435ead5df1d152d8c..c489819ef02c098e11cfb64c05b29587de2b4564 100644
--- a/dumux/freeflow/stokes/stokesproperties.hh
+++ b/dumux/freeflow/stokes/stokesproperties.hh
@@ -27,7 +27,7 @@
  *
  * \file
  *
- * \brief Defines the properties required for the Stokes BOX model.
+ * \brief Defines the properties required for the Stokes box model.
  */
 
 #ifndef DUMUX_STOKESPROPERTIES_HH
@@ -61,7 +61,7 @@ NEW_PROP_TAG(FluidState);
 NEW_PROP_TAG(StabilizationAlpha); //!< The parameter for the stabilization
 NEW_PROP_TAG(StabilizationBeta); //!< The parameter for the stabilization at boundaries
 
-NEW_PROP_TAG(PhaseIndex);
+NEW_PROP_TAG(PhaseIndex); //!< A phase index in case that a two-phase fluidsystem is used
 NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
 NEW_PROP_TAG(Scaling); //!Defines Scaling of the model
 }
diff --git a/dumux/freeflow/stokes/stokespropertydefaults.hh b/dumux/freeflow/stokes/stokespropertydefaults.hh
index eca320a22d9df52593d54f5a9dbe83440ad51409..fb0703d4e521f343ee9edd87a605e6cdc83dc277 100644
--- a/dumux/freeflow/stokes/stokespropertydefaults.hh
+++ b/dumux/freeflow/stokes/stokespropertydefaults.hh
@@ -27,7 +27,10 @@
  *
  * \file
  *
- * \brief Defines the properties required for the Stokes box model.
+ * \brief Defines default properties for the Stokes box model.
+ *
+ * These can be overwritten at a different place or
+ * may be replaced by values of the input file.
  */
 
 #ifndef DUMUX_STOKES_PROPERTY_DEFAULTS_HH
@@ -76,7 +79,7 @@ SET_TYPE_PROP(BoxStokes,
               LocalResidual,
               StokesLocalResidual<TypeTag>);
 
-// Use the Stokes specific newton controller for the Stokes model
+//! Use the Stokes specific newton controller for the Stokes model
 SET_PROP(BoxStokes, NewtonController)
 {public:
     typedef StokesNewtonController<TypeTag> type;
@@ -95,7 +98,7 @@ SET_TYPE_PROP(BoxStokes, VolumeVariables, StokesVolumeVariables<TypeTag>);
 //! the FluxVariables property
 SET_TYPE_PROP(BoxStokes, FluxVariables, StokesFluxVariables<TypeTag>);
 
-//! the upwind factor for the mobility.
+//! the upwind factor.
 SET_SCALAR_PROP(BoxStokes, MassUpwindWeight, 1.0);
 
 //! The fluid system to use by default
@@ -107,6 +110,7 @@ public:
     typedef FluidSystems::OneP<Scalar, Fluid> type;
 };
 
+//! The fluid that is used in the single-phase fluidsystem.
 SET_PROP(BoxStokes, Fluid)
 { private:
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -116,7 +120,7 @@ public:
 
 SET_TYPE_PROP(BoxStokes, StokesIndices, StokesCommonIndices<TypeTag>);
 
-//! Choose the type of the employed fluid state
+//! Choose the type of the employed fluid state.
 SET_PROP(BoxStokes, FluidState)
 {
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -126,6 +130,7 @@ public:
 
 };
 
+//! Set the phaseIndex per default to zero (important for two-phase fluidsystems).
 SET_INT_PROP(BoxStokes, PhaseIndex, 0);
 
 //
diff --git a/dumux/freeflow/stokes/stokesvolumevariables.hh b/dumux/freeflow/stokes/stokesvolumevariables.hh
index 8fa8012e559e6c1b907013899c66f45dacf1a049..20b9e30eaa7d21eeb24bc176208000ca6a7f1575 100644
--- a/dumux/freeflow/stokes/stokesvolumevariables.hh
+++ b/dumux/freeflow/stokes/stokesvolumevariables.hh
@@ -142,38 +142,34 @@ public:
     { return fluidState_; }
 
     /*!
-     * \brief Returns the mass density of a given phase within the
-     *        control volume.
+     * \brief Returns the mass density \f$\mathrm{[kg/m^3]}\f$ of the fluid within the
+     *        sub-control volume.
      */
     Scalar density() const
     { return fluidState_.density(phaseIdx); }
 
     /*!
-     * \brief Returns the effective pressure of a given phase within
-     *        the control volume.
+     * \brief Returns the fluid pressure \f$\mathrm{[Pa]}\f$ within
+     *        the sub-control volume.
      */
     Scalar pressure() const
     { return fluidState_.pressure(phaseIdx); }
 
     /*!
-     * \brief Returns temperature inside the sub-control volume.
-     *
-     * Note that we assume thermodynamic equilibrium, i.e. the
-     * temperature of the rock matrix and of all fluid phases are
-     * identical.
+     * \brief Returns temperature\f$\mathrm{[T]}\f$ inside the sub-control volume.
      */
     Scalar temperature() const
     { return fluidState_.temperature(phaseIdx); }
 
     /*!
-     * \brief Returns the viscosity of the fluid in
-     *        the control volume.
+     * \brief Returns the viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid in
+     *        the sub-control volume.
      */
     Scalar viscosity() const
     { return fluidState_.viscosity(phaseIdx); }
 
     /*!
-     * \brief Returns the phase velocity
+     * \brief Returns the velocity vector in the sub-control volume.
      */
     const VelocityVector &velocity() const
     { return velocity_; }