diff --git a/doc/handbook/5_models.tex b/doc/handbook/5_models.tex
index 6da2bedb66e2bf88b332d0e85686f01a162dcdc1..6fc5e5319374193c59173003d57a21056bcbfd78 100644
--- a/doc/handbook/5_models.tex
+++ b/doc/handbook/5_models.tex
@@ -261,5 +261,5 @@ In comparison to solving the equations monolothically, the sequential structure
 different discretization methods for the different equations. The standard method
 used in the sequential algorithm is a cell-centered finite volume method. Further schemes,
 so far only available for the two-phase pressure equation, are cell-centered finite
-volumes with multi-point flux approximation (MPFA O-method) and mimetic finite differences.
+volumes with multi-point flux approximation (Mpfa-O method) and mimetic finite differences.
 An $h$-adaptive implementation of both sequential algorithms is provided for two dimensions.
diff --git a/doc/handbook/5_spatialdiscretizations.tex b/doc/handbook/5_spatialdiscretizations.tex
index 239ea6723b5cdb7a3205f3cada5cc37b9fe2867a..85ed6d7932e6d7530c18f3f77a12628836f330aa 100644
--- a/doc/handbook/5_spatialdiscretizations.tex
+++ b/doc/handbook/5_spatialdiscretizations.tex
@@ -197,15 +197,15 @@ For the derivation of the finite-volume formulation we integrate the first equat
     \label{eq:ellipticIntegrated}
 \end{equation}
 
-Splitting the control volume boundary $\partial K$ into a finite number of faces $\sigma \subset \partial K$ (such that $\sigma = \overline{K} \cap \overline{L}$ for some neighboring control volume $L$) and replacing the exact fluxes by an approximation, i.e. $F_{K, \sigma} \approx \int_{\sigma} \left( - \mathbf{\Lambda} \nabla u \right) \cdot \mathbf{n} \mathrm{d} \Gamma$, yield
+Splitting the control volume boundary $\partial K$ into a finite number of faces $\sigma \subset \partial K$ (such that $\sigma = \overline{K} \cap \overline{L}$ for some neighboring control volume $L$) and replacing the exact fluxes by an approximation, i.e. $F_{K, \sigma} \approx \int_{\sigma} \left( - \mathbf{\Lambda}_K \nabla u \right) \cdot \mathbf{n} \mathrm{d} \Gamma$ (here $\mathbf{\Lambda}_K$ is the value of $\mathbf{\Lambda}$ associated with control volume $K$), yield
 \begin{equation}
     \sum_{\sigma \subset \partial K} F_{K, \sigma} = Q_K, \quad \forall \, {K \in \mathcal{M}},
 \label{eq:ccdisc}
 \end{equation}
 where $F_{K, \sigma}$ is the discrete flux through face $\sigma$ flowing out of cell $K$ and $Q_K := \int_K q \, \mathrm{d}x$ is the integrated source/sink term. Equation \eqref{eq:ccdisc} is the typical cell-centered finite-volume formulation. 
 Finite-volume schemes differ in the way how the term 
-$(\mathbf{\Lambda} \nabla u ) \cdot \mathbf{n} $ is approximated (i.e. the choice of the fluxes $F_{K, \sigma}$). Using the symmetry of the tensor $\mathbf{\Lambda}$, this term can be rewritten as 
-$\nabla u  \cdot \mathbf{\Lambda}\mathbf{n}$, which corresponds to the directional derivative of $u$ in co-normal direction $\mathbf{\Lambda}\mathbf{n}$. 
+$(\mathbf{\Lambda}_K \nabla u ) \cdot \mathbf{n} $ is approximated (i.e. the choice of the fluxes $F_{K, \sigma}$). Using the symmetry of the tensor $\mathbf{\Lambda}_K$, this term can be rewritten as 
+$\nabla u  \cdot \mathbf{\Lambda}_K\mathbf{n}$, which corresponds to the directional derivative of $u$ in co-normal direction $\mathbf{\Lambda}_K\mathbf{n}$. 
 In the following, the main ideas of the two-point flux approximation and the multi-point flux approximation methods are briefly described. Hereby, we restrict the discussion to the two-dimensional case.
 
 Please also note that other types of equations, e.g. instationary parabolic problems, can be discretized by applying some time discretization scheme to the time derivatives and by using the finite-volume scheme for the flux discretization. For simplicity the discussion is restricted to the elliptic problem \eqref{eq:elliptic}.
@@ -217,7 +217,7 @@ This scheme can be derived by using the conormal decomposition, which reads
 \mathbf{\Lambda}_K \mathbf{n}_{K, \sigma} = t_{K,\sigma} \mathbf{d}_{K,\sigma} + \mathbf{d}^{\bot}_{K,\sigma}, \quad  t_{K,\sigma} = \frac{\mathbf{n}_{K, \sigma}^T \mathbf{\Lambda}_K \mathbf{d}_{K,\sigma} }{\mathbf{d}_{K,\sigma}^T \mathbf{d}_{K,\sigma}}, \; \mathbf{d}^{\bot}_{K,\sigma} = \mathbf{\Lambda}_K \mathbf{n}_{K, \sigma} - t_{K,\sigma} \mathbf{d}_{K,\sigma},
 \label{eq:conormalDecTpfa}
 \end{equation}
-with the distance vector $\mathbf{d}_{K,\sigma} := \mathbf{x}_\sigma - \mathbf{x}_K$ and $\mathbf{d}_{K,\sigma}^T \mathbf{d}^{\bot}_{K,\sigma} = 0$, see Figure \ref{pc:cctpfa} for the used notations. The same can be done for the conormal $\mathbf{\Lambda}_L \mathbf{n}_{L, \sigma}$. The $t_{K,\sigma}$ and $t_{L,\sigma}$ are the transmissibilities associated with the face $\sigma$. These transmissibilities are calculated in \Dumux by using the function \texttt{computeTpfaTransmissibility}.
+with the tensor $\mathbf{\Lambda}_K$ associated with control volume $K$, the distance vector $\mathbf{d}_{K,\sigma} := \mathbf{x}_\sigma - \mathbf{x}_K$ and $\mathbf{d}_{K,\sigma}^T \mathbf{d}^{\bot}_{K,\sigma} = 0$, see Figure \ref{pc:cctpfa} for the used notations. The same can be done for the conormal $\mathbf{\Lambda}_L \mathbf{n}_{L, \sigma}$. The $t_{K,\sigma}$ and $t_{L,\sigma}$ are the transmissibilities associated with the face $\sigma$. These transmissibilities are calculated in \Dumux by using the function \texttt{computeTpfaTransmissibility}.
 
 \begin{figure} [ht]
 \centering
@@ -255,7 +255,7 @@ with $\tau_{K,\sigma} := \mathbf{n}_{K, \sigma} \mathbf{\Lambda}_K\mathbf{n}_{K,
  
 
 \subsubsection{Mpfa Method}\label{cc_mpfa}
-Expressions for the face fluxes $F_{K, \sigma}$ are usually obtained by introducing intermediate face unknowns $u_\sigma$ in addition to the cell unknowns $u_K$ and enforcing the physically motivated continuity of fluxes and continuity of the solution across the faces. For a face $\sigma$ between the two polygons $K$ and $L$ these conditions read:
+Expressions for the face fluxes $F_{K, \sigma}$ are obtained by introducing intermediate face unknowns $u_\sigma$ in addition to the cell unknowns $u_K$ and enforcing the physically motivated continuity of fluxes and continuity of the solution across the faces. For a face $\sigma$ between the two polygons $K$ and $L$ these conditions read:
 \begin{equation}
     \begin{aligned}
         &F_{K, \sigma} + F_{L, \sigma} = 0 \\
@@ -263,7 +263,7 @@ Expressions for the face fluxes $F_{K, \sigma}$ are usually obtained by introduc
         \label{eq:sigmaConditions}
     \end{aligned}
 \end{equation}
-Using these conditions the intermediate face unknowns ${u}_\sigma$ can be eliminated and the fluxes are expressed as a function of the cell unknowns $u_N$ and associated transmissibilities $t^N_{K,\sigma}$:
+Using these conditions, the intermediate face unknowns ${u}_\sigma$ can be eliminated and the fluxes are expressed as a function of the cell unknowns $u_N$ and associated transmissibilities $t^N_{K,\sigma}$:
 
 \begin{equation}
     F_{K,\sigma} = \sum_{N \in \mathcal{S}_{K,\sigma}} t^N_{K,\sigma} u_{N}.
@@ -283,9 +283,9 @@ The main difference between the various finite-volume schemes available is the a
 \end{equation*}
 with $t_{K,\sigma},t_{L,\sigma}$ as defined in equation \eqref{eq:conormalDecTpfa}.
 
-In the following, a multi-point flux approximation method (Mpfa-O method), which was first introduced in \citet{Aavatsmark2002}, is presented. The main difference to the Tpfa scheme is the fact that a consistent discrete gradient is constructed, i.e. the term $\nabla u \cdot \mathbf{d}^{\bot}_{K,\sigma}$ is not neglected.
+In the following, a multi-point flux approximation method (Mpfa-O method), which was introduced in \citet{A3:aavatsmark:2002}, is presented. The main difference to the Tpfa scheme is the fact that a consistent discrete gradient is constructed, i.e. the term $\nabla u \cdot \mathbf{d}^{\bot}_{K,\sigma}$ is not neglected.
 
-For this scheme, a dual grid is created by connecting the barycenters of the cells with the barycenters of the faces ($d=2$) or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control volumes $K^v$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^v$, see Figure \ref{pc:interactionRegion_mpfa}. We allow for piecewise constant $\mathbf{\Lambda}$ (per cell) and construct discrete gradients $\nabla_\mathcal{D}^{K^v} u$ (per sub-control volume $K^v$). 
+For this scheme, a dual grid is created by connecting the barycenters of the cells with the barycenters of the faces ($d=2$) or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control volumes $K^v$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^v$, see Figure \ref{pc:interactionRegion_mpfa}. We allow for piecewise constant $\mathbf{\Lambda}$ (denoted as $\mathbf{\Lambda}_K$ for each cell $K$) and construct discrete gradients $\nabla_\mathcal{D}^{K^v} u$ (per sub-control volume $K^v$). 
 In the following, we restrict our discussion to the two-dimensional setup that is shown in Figure \ref{pc:interactionRegion_mpfa}.  
 Here, the discrete gradients are constructed to be consistent such that the following conditions hold:
 \begin{equation}
@@ -314,7 +314,7 @@ Inserting the discrete gradient, yields
     F_{K, \sigma^v_1} = \omega_{K,\sigma^v_1\sigma^v_1}(u_K - u_{\sigma^v_1}) + \omega_{K,\sigma^v_1 \sigma^v_3}(u_K - u_{\sigma^v_3}),
     \label{eq:discreteFluxRef}
 \end{equation}
-with $(\omega_{K,\sigma^v_1\sigma^v_1},\omega_{K,\sigma^v_1 \sigma^v_3})^T = |\sigma^v_1| \mathbb{D}^{-1}_{K^v}\mathbf{\Lambda}_K \mathbf{n}_{\sigma^v_1}$. 
+with $(\omega_{K,\sigma^v_1\sigma^v_1},\omega_{K,\sigma^v_1 \sigma^v_3})^T = |\sigma^v_1| \mathbb{D}^{-1}_{K^v}\mathbf{\Lambda}_K \mathbf{n}_{\sigma^v_1}$. These values are calculated in \Dumux by using the function \texttt{computeMpfaTransmissibility}.
 \\ \ \\
 To deduce a cell-centered scheme, the introduced face unknowns $u_{\sigma^v_i}$ have to be eliminated. This is done by enforcing flux continuity for each sub-control volume face, i.e.
 \begin{align}
@@ -331,7 +331,6 @@ Inserting these face unknowns into the flux expression \eqref{eq:discreteFluxRef
     \label{eq:FVFluxExpressionSubFace}
 \end{equation}
 for each cell $K$ and sub-control volume face $\sigma^v_i$. 
-In \Dumux the transmissibility vector $\mathbf{t}_{K,\sigma^v_i}$ is returned by the function \texttt{advectionTijSecondaryIv()} or \texttt{advectionTijPrimaryIv()}, depending on the chosen interaction volume type (the primary interaction volume is used by default).
 
 % \subsubsection{NLTPFA}\label{cc_nltpfa}
 % TODO
diff --git a/doc/handbook/CMakeLists.txt b/doc/handbook/CMakeLists.txt
index 01cf2ee849d268170050adcb9ce85ef115e4b8d2..5b3a54ea25982e5640bbc939e4e6cf579352b358 100644
--- a/doc/handbook/CMakeLists.txt
+++ b/doc/handbook/CMakeLists.txt
@@ -23,13 +23,14 @@ set(TEX_INPUTS
 
 set(TEX_IMAGES
   png/box_disc.png
-  png/cc_disc.png
+  png/cctpfa.png
   png/dumux_strucutre_flowchart_horizontal_explained.png
   png/dunedesign.png
   ../logo/dumux_logo_hires_whitebg.png
   png/dalton1.png
   png/dalton2.png
-  pdf/staggered_grid.pdf)
+  pdf/staggered_grid.pdf
+  png/mpfa_iv.png)
 
 dune_add_latex_document(0_dumux-handbook.tex
   BIBFILES dumux-handbook.bib