diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
index 2720ff510c3b3a33e0f0a1357ec0d2107fcc13c8..baeb3b6b87b4ebddc0add67af32093768e865293 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
@@ -35,56 +35,69 @@ namespace Dumux
 //! \ingroup FVPressure2p
 /*!  \brief Finite Volume discretization of a two-phase flow pressure equation of the sequential IMPES model.
  *
- * This model solves equations of the form
- * \f[
- * \phi \left( \rho_w  \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) + \text{div}\, \boldsymbol{v}_{total} = q.
- * \f]
- * The definition of the total velocity \f$\boldsymbol{v}_{total}\f$ depends on the choice of the primary pressure variable.
- * Further, fluids can be assumed to be compressible or incompressible (Property: <tt>EnableCompressibility</tt>).
- * In the incompressible case a wetting \f$(w) \f$ phase pressure as primary variable leads to
- *
- * \f[
- * - \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,
- * \f]
+ * This model implements two-phase flow of two immiscible fluids \f$\alpha \in \{ w, n \}\f$ using
+ * a standard multiphase Darcy approach as the equation for the conservation of momentum, i.e.
+ \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
+ \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right).
+ \f]
+ * By inserting this into the equation for the conservation of the phase mass, one gets
+  \f[
+ \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
+ -
+ \text{div} \left\{
+ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ \right\} = q_\alpha \;.
+ \f]
+
+ * In the incompressible case the phase densities \f$ \varrho_\alpha \f$ can be eliminated from the equations. The two equations
+ * are then added up, and because \f$S_w + S_n = 1\f$, the first term cancels out. This leads to the so-called pressure
+ * equation. For the wetting (\f$w\f$)
+ * phase pressure as primary variable this yields
+ \f[
+ - \text{div}\,  \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_w + f_n \textbf{grad}\,
+   p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q,
+ \f]
  *
- * a non-wetting (\f$ n \f$) phase pressure yields
- * \f[
- *  - \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,
- *  \f]
- * and a global pressure leads to
- * \f[
- * - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\,
- *   p_{global} + \sum f_\alpha \rho_\alpha \, g \, \textbf{grad}\, z\right)\right] = q.
- * \f]
+ * for the non-wetting (\f$ n \f$) phase pressure as primary variable it yields
+ \f[
+  - \text{div}\,  \left[\lambda \boldsymbol K  \left(\textbf{grad}\, p_n - f_w \textbf{grad}\,
+    p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q,
+  \f]
+ * and for the global pressure as primary variable it leads to
+ \f[
+ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\,
+   p_{global} - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q.
+ \f]
  * Here, \f$ p_\alpha \f$ is a phase pressure, \f$ p_ {global} \f$ 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.),
- * \f$ p_c = p_n - p_w \f$ is the capillary pressure, \f$ \boldsymbol K \f$ the absolute permeability,
+ * \f$ p_c = p_n - p_w \f$ is the capillary pressure, \f$ \boldsymbol K \f$ the absolute permeability tensor,
  * \f$ \lambda = \lambda_w +  \lambda_n \f$ the total mobility depending on the
  * saturation (\f$ \lambda_\alpha = k_{r_\alpha} / \mu_\alpha \f$),
  * \f$ f_\alpha = \lambda_\alpha / \lambda \f$ the fractional flow function of a phase,
- * \f$ \rho_\alpha \f$ a phase density, \f$ g \f$ the gravity constant and \f$ q \f$ the source term.
+ * \f$ \varrho_\alpha \f$ a phase density, \f$ {\textbf g} \f$ the gravitational acceleration vector and \f$ q = q_w + q_n \f$ the total source term.
+ * Depending on the primary variable chosen, one of the pressure equations above is solved sequentially together with the conservation
+ * of the phase mass for one phase.
  *
  * For all cases, \f$ p = p_D \f$ on \f$ \Gamma_{Dirichlet} \f$, and \f$ \boldsymbol v_{total} \cdot  \boldsymbol n  = q_N \f$
  * on \f$ \Gamma_{Neumann} \f$.
  *
  * The slightly compressible case is only implemented for phase pressures! In this case for a wetting
- * \f$(w) \f$ phase pressure as primary variable the equations are formulated as
+ * \f$(w) \f$ phase pressure as primary variable the pressure equation is formulated as
  * \f[
- * \phi \left( \rho_w  \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \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,
+ * \phi \left( \varrho_w  \frac{\partial S_w}{\partial t} + \varrho_n \frac{\partial S_n}{\partial t}\right) - \text{div}\,
+ * \left[\lambda \boldsymbol{K} \left(\textbf{grad}\, p_w + f_n \, \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha
+ * {\textbf g}\right)\right] = q,
  * \f]
- * and for a non-wetting (\f$ n \f$) phase pressure as
+ * and for a non-wetting (\f$ n \f$) phase pressure as primary variable as
  *  \f[
- *  \phi \left( \rho_w  \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \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,
+ *  \phi \left( \varrho_w  \frac{\partial S_w}{\partial t} + \varrho_n \frac{\partial S_n}{\partial t}\right) - \text{div}\,
+ * \left[\lambda \boldsymbol{K}  \left(\textbf{grad}\, p_n - f_w \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha
+ * {\textbf g}\right)\right] = q.
  *  \f]
  * In this slightly compressible case the following definitions are valid:
- * \f$ \lambda = \rho_w \lambda_w + \rho_n \lambda_n \f$, \f$ f_\alpha = (\rho_\alpha \lambda_\alpha) / \lambda \f$
+ * \f$ \lambda = \varrho_w \lambda_w + \varrho_n \lambda_n \f$, \f$ f_\alpha = (\varrho_\alpha \lambda_\alpha) / \lambda \f$
  * This model assumes that temporal changes in density are very small and thus terms of temporal derivatives are negligible in the pressure equation.
  * Depending on the formulation the terms including time derivatives of saturations are simplified by inserting  \f$ S_w + S_n = 1 \f$.
  *
@@ -613,7 +626,7 @@ void FVPressure2P<TypeTag>::getSource(EntryType& entry, const Element& element
  *
  * If compressibility is enabled this functions calculates the term
  * \f[
- *      \phi \sum_\alpha \rho_\alpha \frac{\partial S_\alpha}{\partial t} V
+ *      \phi \sum_\alpha \varrho_\alpha \frac{\partial S_\alpha}{\partial t} V
  * \f]
  *
  * In the incompressible case an volume correction term is calculated which corrects
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
index 5a66c5516ea5c58bf5aee4a68c9df79ca1250faf..ccebaafa6f418dba6ce5081922733d429599e0a3 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
@@ -35,12 +35,13 @@ namespace Dumux
 /*! Calculates phase velocities or total velocity from a known pressure field applying a finite volume discretization.
  * The wetting or the non-wetting phase pressure, or the global pressure has to be given as piecewise constant cell values.
  * The phase velocities are calculated following  Darcy's law as
- * \f[
- * \boldsymbol v_\alpha = \lambda_\alpha \boldsymbol K \left(\textbf{grad}\, p_\alpha + \rho_\alpha g  \textbf{grad}\, z \right),
- * \f]
- * where \f$ p_\alpha \f$ denotes the pressure of phase \f$ _\alpha \f$ (wetting or non-wetting),
- * \f$ \boldsymbol K \f$ the absolute permeability, \f$ \lambda_\alpha \f$ the phase mobility,
- * \f$ \rho_\alpha \f$ the phase density and \f$ g \f$ the gravity constant.
+ \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
+ \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right),
+ \f]
+ * where \f$ p_\alpha \f$ denotes the pressure of phase \f$\alpha \in \{ w, n \}\f$,
+ * \f$ \boldsymbol K \f$ the absolute permeability tensor, \f$ \lambda_\alpha \f$ the phase mobility,
+ * \f$ \varrho_\alpha \f$ the phase density and \f$ {\textbf g} \f$ the gravitational acceleration vector.
  * The total velocity is either calculated as sum of the phase velocities
  * \f[
  * \boldsymbol v_{total} = \boldsymbol v_{wetting}+\boldsymbol v_{non-wetting},
@@ -48,7 +49,7 @@ namespace Dumux
  * or with a given global pressure
  * \f[
  * \boldsymbol v_{total} = \lambda_{total} \boldsymbol K \left(\textbf{grad}\,
- *  p_{global} + \sum f_\alpha \rho_\alpha g  \textbf{grad}\, z\right).
+ *  p_{global} - \sum f_\alpha \varrho_\alpha {\textbf g}\right).
  * \f]
  *
  * \tparam TypeTag The Type Tag
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
index 5a9453d02427f5502e9138e6509306dabe8f88ef..1563327874a8dd490361e59f093d619ebb808b44 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
@@ -34,31 +34,31 @@ namespace Dumux
 /*! This model solves equations of the form
  *
  *  \f[
- *  \phi \frac{\partial (\rho_\alpha S_\alpha)}{\partial t} + \text{div}\, (\rho_\alpha \boldsymbol{v_\alpha}) = q_\alpha,
+ *  \phi \frac{\partial (\varrho_\alpha S_\alpha)}{\partial t} + \text{div}\, (\varrho_\alpha \boldsymbol{v_\alpha}) = q_\alpha,
  *  \f]
  *
- *  where \f$ S_\alpha \f$ is the saturation of phase \f$ \alpha \f$ (wetting \f$(w) \f$,
- *  non-wetting \f$(n) \f$) and \f$ \boldsymbol v_\alpha \f$ is the phase velocity defined by
+ *  where \f$ S_\alpha \f$ is the saturation of phase \f$\alpha \in \{ w, n \}\f$
+ *  and \f$ \boldsymbol v_\alpha \f$ is the phase velocity defined by
  *  the multi-phase Darcy equation.
  *  If a phase velocity is reconstructed from the pressure solution it can be directly inserted into
  *  the previous equation. In the incompressible case the equation is further divided by the phase density
- *  \f$ \rho_\alpha \f$. If a total velocity is reconstructed the saturation equation is reformulated into:
+ *  \f$ \varrho_\alpha \f$. If a total velocity is reconstructed the saturation equation is reformulated into:
  *
  * \f[
  *  \phi \frac{\partial S_w}{\partial t} + f_w \text{div}\, \boldsymbol{v}_{t} + f_w \lambda_n \boldsymbol{K}\left(\textbf{grad}\,
- *  p_c + (\rho_n-\rho_w) \, g \, \textbf{grad} z \right)= q_\alpha,
+ *  p_c - (\varrho_n-\varrho_w) {\textbf g} \right)= q_\alpha,
  * \f]
  * to get a wetting phase saturation or
  * \f[
  * \phi \frac{\partial S_n}{\partial t} + f_n \text{div}\, \boldsymbol{v}_{t} - f_n \lambda_w \boldsymbol{K}\left(\textbf{grad}\,
- * p_c + (\rho_n-\rho_w) \, g \, \textbf{grad} z \right)= q_\alpha,
+ * p_c - (\varrho_n-\varrho_w) {\textbf g} \right)= q_\alpha,
  * \f]
  * if the non-wetting phase saturation is the primary transport variable.
  *
  *  The total velocity formulation is only implemented for incompressible fluids and \f$ f_\alpha \f$
  *  is the fractional flow function, \f$ \lambda_\alpha \f$ is the mobility, \f$ \boldsymbol K \f$
- *  the absolute permeability,\f$ p_c \f$ the capillary pressure, \f$ \rho \f$ the fluid density,
- *  \f$ g \f$ the gravity constant, and \f$ q \f$ the source term.
+ *  the absolute permeability tensor,\f$ p_c \f$ the capillary pressure, \f$ \varrho_\alpha \f$ the phase density,
+ *  \f$ {\textbf g} \f$ the gravitational acceleration vector, and \f$ q_\alpha \f$ the source term.
  *
  *
  *  In the IMPES models the default setting is: