diff --git a/doc/handbook/ModelDescriptions/2p2cboxmodel.tex b/doc/handbook/ModelDescriptions/2p2cboxmodel.tex
index 6e39c9d747e45209b70b8a1d2eebe8c7f7ab79d7..23d1ad2f8f31b3ac548b1fc9a43495ac8956e4bb 100644
--- a/doc/handbook/ModelDescriptions/2p2cboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2p2cboxmodel.tex
@@ -6,14 +6,14 @@
 
 \-This model implements two-\/phase two-\/component flow of two compressible and partially miscible fluids $\alpha \in \{ w, n \}$ composed of the two components $\kappa \in \{ w, a \}$. \-The standard multiphase \-Darcy approach is used as the equation for the conservation of momentum\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \]
 
-\-By inserting this into the equations for the conservation of the components, one gets one transport equation for each component \begin{eqnarray} && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )} {\partial t} - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} (\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ &-& \sum_\alpha \text{div} \left\{{\bf D}_{\alpha, pm}^\kappa \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\} - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, , \alpha \in \{w, g\} \end{eqnarray}
+\-By inserting this into the equations for the conservation of the components, one gets one transport equation for each component \begin{eqnarray*} && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )} {\partial t} - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} (\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ &-& \sum_\alpha \text{div} \left\{{\bf D}_{\alpha, pm}^\kappa \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\} - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, , \alpha \in \{w, g\} \end{eqnarray*}
 
 \-This is discretized using a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit \-Euler method as temporal discretization.
 
 \-By using constitutive relations for the capillary pressure $p_c = p_n - p_w$ and relative permeability $k_{r\alpha}$ and taking advantage of the fact that $S_w + S_n = 1$ and $X^\kappa_w + X^\kappa_n = 1$, the number of unknowns can be reduced to two. \-The used primary variables are, like in the two-\/phase model, either $p_w$ and $S_n$ or $p_n$ and $S_w$. \-The formulation which ought to be used can be specified by setting the {\ttfamily \-Formulation} property to either \-Two\-P\-Two\-C\-Indices\-::p\-Ws\-N or \-Two\-P\-Two\-C\-Indices\-::p\-Ns\-W. \-By default, the model uses $p_w$ and $S_n$. \-Moreover, the second primary variable depends on the phase state, since a primary variable switch is included. \-The phase state is stored for all nodes of the system. \-Following cases can be distinguished\-:
 \begin{itemize}
 \item \-Both phases are present\-: \-The saturation is used (either $S_n$ or $S_w$, dependent on the chosen {\ttfamily \-Formulation}), as long as $ 0 < S_\alpha < 1$.
-\item \-Only wetting phase is present\-: \-The mass fraction of, e.\-g., air in the wetting phase $X^a_w$ is used, as long as the maximum mass fraction is not exceeded ( $X^a_w<X^a_{w,max}$)
-\item \-Only non-\/wetting phase is present\-: \-The mass fraction of, e.\-g., water in the non-\/wetting phase, $X^w_n$, is used, as long as the maximum mass fraction is not exceeded ( $X^w_n<X^w_{n,max}$)
+\item \-Only wetting phase is present\-: \-The mass fraction of, e.\-g., air in the wetting phase $X^a_w$ is used, as long as the maximum mass fraction is not exceeded $(X^a_w<X^a_{w,max})$
+\item \-Only non-\/wetting phase is present\-: \-The mass fraction of, e.\-g., water in the non-\/wetting phase, $X^w_n$, is used, as long as the maximum mass fraction is not exceeded $(X^w_n<X^w_{n,max})$
 \end{itemize}
 
diff --git a/doc/handbook/ModelDescriptions/2p2cniboxmodel.tex b/doc/handbook/ModelDescriptions/2p2cniboxmodel.tex
index c77d992b8a16b6e6d8a52b5aea6feb401af64d56..fa612e102438bdceecc4e23f6cd35dbdfa908ee0 100644
--- a/doc/handbook/ModelDescriptions/2p2cniboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2p2cniboxmodel.tex
@@ -8,7 +8,7 @@
 
 \-This is discretized using a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit \-Euler method as temporal discretization.
 
-\-By using constitutive relations for the capillary pressure $p_c = p_n - p_w$ and relative permeability $k_{r\alpha}$ and taking advantage of the fact that $S_w + S_n = 1$ and $X^\kappa_w + X^\kappa_n = 1$, the number of unknowns can be reduced to two. \-If both phases are present the primary variables are, like in the nonisothermal two-\/phase model, either $p_w$, $S_n$ and temperature or $p_n$, $S_w$ and temperature. \-The formulation which ought to be used can be specified by setting the {\ttfamily \-Formulation} property to either {\ttfamily \-Two\-P\-Two\-Indices\-::p\-Ws\-N} or {\ttfamily \-Two\-P\-Two\-C\-Indices\-::p\-Ns\-W}. \-By default, the model uses $p_w$ and $S_n$. \-In case that only one phase (nonwetting or wetting phase) is present the second primary variable represents a mass fraction. \-The correct assignment of the second primary variable is performed by a phase state dependent primary variable switch. \-The phase state is stored for all nodes of the system. \-The following cases can be distinguished\-:
+\-By using constitutive relations for the capillary pressure $p_c = p_n - p_w$ and relative permeability $k_{r\alpha}$ and taking advantage of the fact that $S_w + S_n = 1$ and $X^\kappa_w + X^\kappa_n = 1$, the number of unknowns can be reduced to two. \-If both phases are present the primary variables are, like in the nonisothermal two-\/phase model, either $p_w$, $S_n$ and temperature or $p_n$, $S_w$ and temperature. \-The formulation which ought to be used can be specified by setting the {\ttfamily \-Formulation} property to either {\ttfamily \-Two\-P\-Two\-C\-Indices\-::p\-Ws\-N} or {\ttfamily \-Two\-P\-Two\-C\-Indices\-::p\-Ns\-W}. \-By default, the model uses $p_w$ and $S_n$. \-In case that only one phase (nonwetting or wetting phase) is present the second primary variable represents a mass fraction. \-The correct assignment of the second primary variable is performed by a phase state dependent primary variable switch. \-The phase state is stored for all nodes of the system. \-The following cases can be distinguished\-:
 \begin{itemize}
 \item \-Both phases are present\-: \-The saturation is used (either $S_n$ or $S_w$, dependent on the chosen formulation).
 \item \-Only wetting phase is present\-: \-The mass fraction of air in the wetting phase $X^a_w$ is used.
diff --git a/doc/handbook/ModelDescriptions/2pboxmodel.tex b/doc/handbook/ModelDescriptions/2pboxmodel.tex
index 6817d2a1fbf9873f93ece43357105ced0bb0b1a0..ad175aea3764f79aa3067e5aa5c763f33cfee890 100644
--- a/doc/handbook/ModelDescriptions/2pboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pboxmodel.tex
@@ -8,7 +8,7 @@
 
 \-By inserting this into the equation for the conservation of the phase mass, one gets \[ \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t} - \text{div} \left\{ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \right\} - q_\alpha = 0 \;, \]
 
-\-This equations are by a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit \-Euler method as time discretization.
+\-These equations are discretized by a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit \-Euler method as time discretization.
 
 \-By using constitutive relations for the capillary pressure $p_c = p_n - p_w$ and relative permeability $k_{r\alpha}$ and taking advantage of the fact that $S_w + S_n = 1$, the number of unknowns can be reduced to two. \-Currently the model supports choosing either $p_w$ and $S_n$ or $p_n$ and $S_w$ as primary variables. \-The formulation which ought to be used can be specified by setting the {\ttfamily \-Formulation} property to either {\ttfamily \-Two\-P\-Common\-Indices\-::p\-Ws\-N} or {\ttfamily \-Two\-P\-Common\-Indices\-::p\-Ns\-W}. \-By default, the model uses $p_w$ and $S_n$.
 
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
index 705badf061335d0ca398713f098ae6cb6a14631b..42c8d6d8c8c11bb1ded93f37b7151d705faec635 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
@@ -4,15 +4,15 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\-This model solves equations of the form \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) + \textbf{div}\, \boldsymbol{v}_{total} = q. \] \-The definition of the total velocity $\boldsymbol{v}_{total}$ depends on the choice of the primary pressure variable. \-Further, fluids can be assumed to be compressible or incompressible (\-Property\-: {\ttfamily \-Enable\-Compressibility}). \-In the incompressible case a wetting ( $ w $) phase pressure as primary variable leads to
+\-This model solves equations of the form \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) + \textbf{div}\, \boldsymbol{v}_{total} = q. \] \-The definition of the total velocity $\boldsymbol{v}_{total}$ depends on the choice of the primary pressure variable. \-Further, fluids can be assumed to be compressible or incompressible (\-Property\-: {\ttfamily \-Enable\-Compressibility}). \-In the incompressible case a wetting $(w) $ phase pressure as primary variable leads to
 
 \[ - \textbf{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 \[ - \textbf{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 \[ - \textbf{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 \[ - \textbf{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 \[ - \textbf{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.
 
 \-For all cases, $ p = p_D $ on $ \Gamma_{Dirichlet} $, and $ \boldsymbol v_{total} \cdot \boldsymbol n = q_N $ on $ \Gamma_{Neumann} $.
 
-\-The slightly compressible case is only implemented for phase pressures! \-In this case for a wetting ( $ w $) phase pressure as primary variable the equations are formulated as \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \textbf{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, \] and for a non-\/wetting ( $ n $) phase pressure as \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \textbf{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, \] \-In this slightly compressible case the following definitions are valid\-: $ \lambda = \rho_w \lambda_w + \rho_n \lambda_n $, $ f_\alpha = (\rho_\alpha \lambda_\alpha) / \lambda $ \-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 $ S_w + S_n = 1 $.
+\-The slightly compressible case is only implemented for phase pressures! \-In this case for a wetting $(w) $ phase pressure as primary variable the equations are formulated as \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \textbf{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, \] and for a non-\/wetting ( $ n $) phase pressure as \[ \phi \left( \rho_w \frac{\partial S_w}{\partial t} + \rho_n \frac{\partial S_n}{\partial t}\right) - \textbf{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, \] \-In this slightly compressible case the following definitions are valid\-: $ \lambda = \rho_w \lambda_w + \rho_n \lambda_n $, $ f_\alpha = (\rho_\alpha \lambda_\alpha) / \lambda $ \-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 $ S_w + S_n = 1 $.
 
 \-In the \-I\-M\-P\-E\-S models the default setting is\-:
 
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
index e794e3f42d402398f19e21bfa81b02758a3e8ccc..62bf3d6ae1145f7a2f673396bb8b2bb0ae9544e7 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
@@ -8,7 +8,7 @@
 
 \[ \phi \frac{\partial (\rho_\alpha S_\alpha)}{\partial t} + \textbf{div}\, (\rho_\alpha \boldsymbol{v_\alpha}) = q_\alpha, \]
 
-where $ S_\alpha $ is the saturation of phase alpha (wetting ( $ w $), non-\/wetting ( $ n $)) and $ \boldsymbol v_\alpha $ 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 in the previous equation. \-In the incompressible case the equation is further divided by the phase density $ \rho_\alpha $. \-If a total velocity is reconstructed the saturation equation is reformulated into\-:
+where $ S_\alpha $ is the saturation of phase $ \alpha $ (wetting $(w) $, non-\/wetting $(n) $) and $ \boldsymbol v_\alpha $ 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 $ \rho_\alpha $. \-If a total velocity is reconstructed the saturation equation is reformulated into\-:
 
 \[ \phi \frac{\partial S_w}{\partial t} + f_w \textbf{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, \] to get a wetting phase saturation or \[ \phi \frac{\partial S_n}{\partial t} + f_n \textbf{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, \] if the non-\/wetting phase saturation is the primary transport variable.
 
diff --git a/doc/handbook/ModelDescriptions/2pniboxmodel.tex b/doc/handbook/ModelDescriptions/2pniboxmodel.tex
index 05bf28e75c625bfabee058e8b0054b6d4aa7d997..9e1995cf51c40257424151621533907e9a9502e3 100644
--- a/doc/handbook/ModelDescriptions/2pniboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pniboxmodel.tex
@@ -8,7 +8,7 @@
 
 \-For the energy balance, local thermal equilibrium is assumed. \-This results in one energy conservation equation for the porous solid matrix and the fluids\-:
 
-\begin{align*} \frac{\partial \phi \sum_alpha \varrho_\alpha u_\alpha S_\alpha}{\partial t} & + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t} - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \textbf{grad}\,p_\alpha - \varrho_\alpha \mbox{\bf g} \right) \right\} \\ & - \text{div} \left(\lambda_{pm} \textbf{grad} \, T \right) - q^h = 0, \qquad \alpha \in \{w, n\} \;, \end{align*} where $h_\alpha$ is the specific enthalpy of a fluid phase $\alpha$ and $u_\alpha = h_\alpha - p_\alpha/\varrho_\alpha$ is the specific internal energy of the phase.
+\begin{align*} \frac{\partial \phi \sum_\alpha \varrho_\alpha u_\alpha S_\alpha}{\partial t} & + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t} - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \textbf{grad}\,p_\alpha - \varrho_\alpha \mbox{\bf g} \right) \right\} \\ & - \text{div} \left(\lambda_{pm} \textbf{grad} \, T \right) - q^h = 0, \qquad \alpha \in \{w, n\} \;, \end{align*} where $h_\alpha$ is the specific enthalpy of a fluid phase $\alpha$ and $u_\alpha = h_\alpha - p_\alpha/\varrho_\alpha$ is the specific internal energy of the phase.
 
 \-The equations are discretized using a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit \-Euler method as time discretization.
 
diff --git a/doc/handbook/ModelDescriptions/3p3cboxmodel.tex b/doc/handbook/ModelDescriptions/3p3cboxmodel.tex
index 5bd6b88d77c698c649281eeb7e95a73f08ce1550..cc8a58ce199b0ecf65359e1e12df6f98694f8406 100644
--- a/doc/handbook/ModelDescriptions/3p3cboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/3p3cboxmodel.tex
@@ -6,7 +6,7 @@
 
 \-This model implements three-\/phase three-\/component flow of three fluid phases $\alpha \in \{ water, gas, NAPL \}$ each composed of up to three components $\kappa \in \{ water, air, contaminant \}$. \-The standard multiphase \-Darcy approach is used as the equation for the conservation of momentum\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \]
 
-\-By inserting this into the equations for the conservation of the components, one transport equation for each component is obtained as \begin{eqnarray} && \phi \frac{\partial (\sum_\alpha \varrho_{\text{mol}, \alpha} x_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_{\text{mol}, \alpha} x_\alpha^\kappa \mbox{\bf K} (\text{grad}\, p_\alpha - \varrho_{\text{mass}, \alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_{pm}^\kappa \varrho_{\text{mol}, \alpha } \text{grad}\, x_\alpha^\kappa \right\} - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha \end{eqnarray}
+\-By inserting this into the equations for the conservation of the components, one transport equation for each component is obtained as \begin{eqnarray*} && \phi \frac{\partial (\sum_\alpha \varrho_{\text{mol}, \alpha} x_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_{\text{mol}, \alpha} x_\alpha^\kappa \mbox{\bf K} (\text{grad}\, p_\alpha - \varrho_{\text{mass}, \alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_{pm}^\kappa \varrho_{\text{mol}, \alpha } \text{grad}\, x_\alpha^\kappa \right\} - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha \end{eqnarray*}
 
 \-Note that these balance equations are molar.
 
@@ -16,11 +16,11 @@
 
 \-The used primary variables are dependent on the locally present fluid phases \-An adaptive primary variable switch is included. \-The phase state is stored for all nodes of the system. \-The following cases can be distinguished\-:
 \begin{itemize}
-\item \-All three phases are present\-: \-Primary variables are two saturations ( $S_w$ and $S_n$, and a pressure, in this case $p_g$.
-\item \-Only the water phase is present\-: \-Primary variables are now the mole fractions of air and contaminant in the water phase ( $x_w^a$ and $x_w^c$), as well as the gas pressure, which is, of course, in a case where only the water phase is present, just the same as the water pressure.
-\item \-Gas and \-N\-A\-P\-L phases are present\-: \-Primary variables ( $S_n$, $x_g^w$, $p_g$).
-\item \-Water and \-N\-A\-P\-L phases are present\-: \-Primary variables ( $S_n$, $x_w^a$, $p_g$).
-\item \-Only gas phase is present\-: \-Primary variables ( $x_g^w$, $x_g^c$, $p_g$).
-\item \-Water and gas phases are present\-: \-Primary variables ( $S_w$, $x_w^g$, $p_g$).
+\item \-All three phases are present\-: \-Primary variables are two saturations $(S_w$ and $S_n)$, and a pressure, in this case $p_g$.
+\item \-Only the water phase is present\-: \-Primary variables are now the mole fractions of air and contaminant in the water phase $(x_w^a$ and $x_w^c)$, as well as the gas pressure, which is, of course, in a case where only the water phase is present, just the same as the water pressure.
+\item \-Gas and \-N\-A\-P\-L phases are present\-: \-Primary variables $(S_n$, $x_g^w$, $p_g)$.
+\item \-Water and \-N\-A\-P\-L phases are present\-: \-Primary variables $(S_n$, $x_w^a$, $p_g)$.
+\item \-Only gas phase is present\-: \-Primary variables $(x_g^w$, $x_g^c$, $p_g)$.
+\item \-Water and gas phases are present\-: \-Primary variables $(S_w$, $x_w^g$, $p_g)$.
 \end{itemize}
 
diff --git a/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex b/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex
index 06cbf16327230451f95c06bb2229ebe75ef4fa2b..92161de806362071ca758e969212fccd8313154b 100644
--- a/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex
@@ -16,11 +16,11 @@
 
 \-The used primary variables are dependent on the locally present fluid phases \-An adaptive primary variable switch is included. \-The phase state is stored for all nodes of the system. \-The following cases can be distinguished\-:
 \begin{itemize}
-\item \-All three phases are present\-: \-Primary variables are two saturations ( $S_w$ and $S_n$, a pressure, in this case $p_g$, and the temperature $T$.
-\item \-Only the water phase is present\-: \-Primary variables are now the mole fractions of air and contaminant in the water phase ( $x_w^a$ and $x_w^c$), as well as temperature and the gas pressure, which is, of course, in a case where only the water phase is present, just the same as the water pressure.
-\item \-Gas and \-N\-A\-P\-L phases are present\-: \-Primary variables ( $S_n$, $x_g^w$, $p_g$, $T$).
-\item \-Water and \-N\-A\-P\-L phases are present\-: \-Primary variables ( $S_n$, $x_w^a$, $p_g$, $T$).
-\item \-Only gas phase is present\-: \-Primary variables ( $x_g^w$, $x_g^c$, $p_g$, $T$).
-\item \-Water and gas phases are present\-: \-Primary variables ( $S_w$, $x_w^g$, $p_g$, $T$).
+\item \-All three phases are present\-: \-Primary variables are two saturations $(S_w$ and $S_n)$, a pressure, in this case $p_g$, and the temperature $T$.
+\item \-Only the water phase is present\-: \-Primary variables are now the mole fractions of air and contaminant in the water phase $(x_w^a$ and $x_w^c)$, as well as temperature and the gas pressure, which is, of course, in a case where only the water phase is present, just the same as the water pressure.
+\item \-Gas and \-N\-A\-P\-L phases are present\-: \-Primary variables $(S_n$, $x_g^w$, $p_g$, $T)$.
+\item \-Water and \-N\-A\-P\-L phases are present\-: \-Primary variables $(S_n$, $x_w^a$, $p_g$, $T)$.
+\item \-Only gas phase is present\-: \-Primary variables $(x_g^w$, $x_g^c$, $p_g$, $T)$.
+\item \-Water and gas phases are present\-: \-Primary variables $(S_w$, $x_w^g$, $p_g$, $T)$.
 \end{itemize}
 
diff --git a/doc/handbook/ModelDescriptions/mpncboxmodel.tex b/doc/handbook/ModelDescriptions/mpncboxmodel.tex
index 4aa4f60dfcfabeac155308d590f11a446490b603..b47e9c9c9f2e93b857e1bd2f43b5494cada1a81a 100644
--- a/doc/handbook/ModelDescriptions/mpncboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/mpncboxmodel.tex
@@ -8,13 +8,13 @@
 
 \-The standard multi-\/phase \-Darcy law is used as the equation for the conservation of momentum\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K} \left( \text{grad}\left(p_\alpha - \varrho_{\alpha} g\right) \right) \]
 
-\-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 \varrho_\alpha x_\alpha^\kappa S_\alpha}{\partial t} - \mathrm{div}\; \left\{ \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K} \mathbf{grad}\left( p_\alpha - \varrho_{\alpha} g\right) \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 \]
+\-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 \varrho_\alpha x_\alpha^\kappa S_\alpha}{\partial t} - \mathrm{div}\; \left\{ \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K} \text{grad}\left( p_\alpha - \varrho_{\alpha} g\right) \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 \]
 
 \-For the missing $M$ model assumptions, the model assumes that if a fluid phase is not present, the sum of the mole fractions of this fluid phase is smaller than $1$, i.\-e. \[ \forall \alpha: S_\alpha = 0 \implies \sum_\kappa x_\alpha^\kappa \leq 1 \]
 
 \-Also, if a fluid phase may be present at a given spatial location its saturation must be positive\-: \[ \forall \alpha: \sum_\kappa x_\alpha^\kappa = 1 \implies S_\alpha \geq 0 \]
 
-\-Since at any given spatial location, a phase is always either present or not present, the one of the strict equalities on the right hand side is always true, i.\-e. \[ \forall \alpha: S_\alpha \left( \sum_\kappa x_\alpha^\kappa - 1 \right) = 0 \] always holds.
+\-Since at any given spatial location, a phase is always either present or not present, one of the strict equalities on the right hand side is always true, i.\-e. \[ \forall \alpha: S_\alpha \left( \sum_\kappa x_\alpha^\kappa - 1 \right) = 0 \] always holds.
 
 \-These three equations constitute a non-\/linear complementarity problem, which can be solved using so-\/called non-\/linear complementarity functions $\Phi(a, b)$ which have the property \[\Phi(a,b) = 0 \iff a \geq0 \land b \geq0 \land a \cdot b = 0 \]