diff --git a/doc/doxygen/modules b/doc/doxygen/modules
index 19f3c81eea9dbd5de0e7ca095277e4da473b088b..5648fd6cf1505f5b1f81405cb737f834ea49930d 100644
--- a/doc/doxygen/modules
+++ b/doc/doxygen/modules
@@ -201,6 +201,24 @@
          *
          * \copydetails Dumux::Stokes2cniModel
          */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup ElasticBoxModel Linear elastic 
+         *
+         * \copydetails Dumux::ElasticModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup ElOnePTwoCBoxModel One-phase two component linear elastic 
+         *
+         * \copydetails Dumux::ElOnePTwoCModel
+         */
+        /*!
+         * \ingroup ImplicitModels
+         * \defgroup ElTwoPBoxModel Two-phase immiscible linear elastic 
+         *
+         * \copydetails Dumux::ElTwoPModel
+         */
     /*!
      * \ingroup ImplicitModel
      * \defgroup ImplicitBaseProblems Base Problems
diff --git a/doc/handbook/ModelDescriptions/1p2cboxmodel.tex b/doc/handbook/ModelDescriptions/1p2cimplicitmodel.tex
similarity index 83%
rename from doc/handbook/ModelDescriptions/1p2cboxmodel.tex
rename to doc/handbook/ModelDescriptions/1p2cimplicitmodel.tex
index 4a1a41b630c4dc07087667ac877b4f78c3ed9dd5..c078d894f45ab7cfcf8aabe56265af6373b5d732 100644
--- a/doc/handbook/ModelDescriptions/1p2cboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/1p2cimplicitmodel.tex
@@ -10,7 +10,7 @@ Gravity can be enabled or disabled via the property system. By inserting this in
 
 The transport of the components $\kappa \in \{ w, a \}$ is described by the following equation\-: \[ \phi \frac{ \partial \varrho X^\kappa}{\partial t} - \text{div} \left\lbrace \varrho X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p - \varrho {\textbf g} \right) + \varrho D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa \right\rbrace = q. \]
 
-All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization.
+All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization. The model is able to use either mole or mass fractions. The property use\-Moles can be set to either true or false in the problem file. Make sure that the according units are used in the problem setup. use\-Moles is set to true by default.
 
 The primary variables are the pressure $p$ and the mole or mass fraction of dissolved component $x$.
 
diff --git a/doc/handbook/ModelDescriptions/1pboxmodel.tex b/doc/handbook/ModelDescriptions/1pimplicitmodel.tex
similarity index 100%
rename from doc/handbook/ModelDescriptions/1pboxmodel.tex
rename to doc/handbook/ModelDescriptions/1pimplicitmodel.tex
diff --git a/doc/handbook/ModelDescriptions/2p2cboxmodel.tex b/doc/handbook/ModelDescriptions/2p2cimplicitmodel.tex
similarity index 56%
rename from doc/handbook/ModelDescriptions/2p2cboxmodel.tex
rename to doc/handbook/ModelDescriptions/2p2cimplicitmodel.tex
index b80ad8c97316e69d17bb1cdbe939328cfadaaa38..41663d06e6ad506ca2d4676cac6ea537f2cafdea 100644
--- a/doc/handbook/ModelDescriptions/2p2cboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2p2cimplicitmodel.tex
@@ -4,16 +4,16 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-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(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \]
+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} \mathbf{K} \left(\textbf{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} (\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ &-& \sum_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_{\alpha} \frac{M^\kappa}{M_\alpha} \textbf{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 \frac{M^\kappa}{M_\alpha} x_\alpha^\kappa S_\alpha )} {\partial t} - \sum_\alpha \text{div} \left\{ \varrho_\alpha \frac{M^\kappa}{M_\alpha} x_\alpha^\kappa \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} (\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ &-& \sum_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_{\alpha} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\} - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, , \alpha \in \{w, g\} \end{eqnarray*}
 
 All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume 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$ 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\-:
+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. The model is able to use either mole or mass fractions. The property use\-Moles can be set to either true or false in the problem file. Make sure that the according units are used in the problem setup. use\-Moles is set to true by default. 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 mole fraction of, e.\-g., air in the wetting phase $x^a_w$ is used, as long as the maximum mole fraction is not exceeded $(x^a_w<x^a_{w,max})$
+\item Only non-\/wetting phase is present\-: The mole fraction of, e.\-g., water in the non-\/wetting phase, $x^w_n$, is used, as long as the maximum mole 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/2p2cniimplicitmodel.tex
similarity index 69%
rename from doc/handbook/ModelDescriptions/2p2cniboxmodel.tex
rename to doc/handbook/ModelDescriptions/2p2cniimplicitmodel.tex
index 43f4241dda539aa96d12b54cda116f2f6fc39ad1..c2a08cc4b0d80c27281c51c2bef90817f69bda60 100644
--- a/doc/handbook/ModelDescriptions/2p2cniboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2p2cniimplicitmodel.tex
@@ -4,7 +4,7 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements a non-\/isothermal two-\/phase flow of two compressible and partly miscible fluids $\alpha \in \{ w, n \}$. Thus each component $\kappa \in \{ w, a \}$ can be present in each phase. Using the standard multiphase Darcy approach a mass balance equation is solved\-: \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} (\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\ &-& \sum_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_{\alpha} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\} - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, , \alpha \in \{w, n\} \end{eqnarray*} For the energy balance, local thermal equilibrium is assumed which results in one energy conservation equation for the porous solid matrix and the fluids\-: \begin{eqnarray*} && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\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 \mathbf{g} \right) \right\} \\ &-& \text{div} \left( \lambda_\text{pm} \textbf{grad} \, T \right) - q^h = 0 \qquad \alpha \in \{w, n\} \end{eqnarray*}
+This model implements a non-\/isothermal two-\/phase flow of two compressible and partly miscible fluids $\alpha \in \{ w, n \}$. Thus each component $\kappa \in \{ w, a \}$ can be present in each phase. Using the standard multiphase Darcy approach a mass balance equation is solved\-: \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} \mathbf{K} (\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\ &-& \sum_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_{\alpha} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\} - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, , \alpha \in \{w, n\} \end{eqnarray*} For the energy balance, local thermal equilibrium is assumed which results in one energy conservation equation for the porous solid matrix and the fluids\-: \begin{eqnarray*} && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\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 \mathbf{g} \right) \right\} \\ &-& \text{div} \left( \lambda_\text{pm} \textbf{grad} \, T \right) - q^h = 0 \qquad \alpha \in \{w, n\} \end{eqnarray*}
 
 All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization.
 
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
index 5a11f6346a7ba91084b7f6a6111b20591baa27c6..5e8ce73d1c353274f344058299c21b424fd14424 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledpressuremodel.tex
@@ -18,11 +18,7 @@ In the I\-M\-P\-E\-S models the default setting is\-:
 
 
 \begin{itemize}
-\item formulation\-: $ p_w-S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00080_a04294fbcf0af5328016a160dbd8bfff9}{Decoupled\-Two\-P\-Common\-Indices\-::pw\-Sw}})
-\end{itemize}
-
-
-\begin{itemize}
+\item formulation\-: $ p_w-S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00096_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
 \item compressibility\-: disabled (Property\-: {\itshape Enable\-Compressibility} set to {\itshape false})
 \end{itemize}
 
diff --git a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
index 31d3e581ca59368d862041c57c0d36b2d5d8f3ca..7f121c5bb21ac986e07d7011c3ea45bb6a31ef87 100644
--- a/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pdecoupledsaturationmodel.tex
@@ -16,7 +16,7 @@ The total velocity formulation is only implemented for incompressible fluids and
 
 In the I\-M\-P\-E\-S models the default setting is\-:
 
-formulation\-: $ p_w $ -\/ $ S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00080_a04294fbcf0af5328016a160dbd8bfff9}{Decoupled\-Two\-P\-Common\-Indices\-::pw\-Sw}})
+formulation\-: $ p_w $ -\/ $ S_w $ (Property\-: {\itshape Formulation} defined as {\itshape \hyperlink{a00096_a601a847774d6e1b2e2a2b469f70c3f22}{Decoupled\-Two\-P\-Common\-Indices\-::pwsw}})
 
 compressibility\-: disabled (Property\-: {\itshape Enable\-Compressibility} set to {\itshape false})
 
diff --git a/doc/handbook/ModelDescriptions/2pdfmimplicitmodel.tex b/doc/handbook/ModelDescriptions/2pdfmimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..eae81f3749d36da23724abe3801c9f0677dba05a
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/2pdfmimplicitmodel.tex
@@ -0,0 +1,14 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements two-\/phase flow of two immiscible fluids $\alpha \in \{ w, n \}$ using a standard multiphase Darcy approach as the equation for the conservation of momentum, i.\-e. \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right) \]
+
+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} \mathbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \right\} - q_\alpha = 0 \;, \]
+
+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/2pboxmodel.tex b/doc/handbook/ModelDescriptions/2pimplicitmodel.tex
similarity index 92%
rename from doc/handbook/ModelDescriptions/2pboxmodel.tex
rename to doc/handbook/ModelDescriptions/2pimplicitmodel.tex
index b7727b6b5156467f809d6ebb6fcc9c56088aebb7..c686b7e69f8815120da1054aacb999157ccc64a4 100644
--- a/doc/handbook/ModelDescriptions/2pboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/2pimplicitmodel.tex
@@ -6,7 +6,7 @@
 
 This model implements two-\/phase flow of two immiscible fluids $\alpha \in \{ w, n \}$ using a standard multiphase Darcy approach as the equation for the conservation of momentum, i.\-e. \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right) \]
 
-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 \;, \]
+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} \mathbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \right\} - q_\alpha = 0 \;, \]
 
 All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization.
 
diff --git a/doc/handbook/ModelDescriptions/2pniboxmodel.tex b/doc/handbook/ModelDescriptions/2pniimplicitmodel.tex
similarity index 100%
rename from doc/handbook/ModelDescriptions/2pniboxmodel.tex
rename to doc/handbook/ModelDescriptions/2pniimplicitmodel.tex
diff --git a/doc/handbook/ModelDescriptions/3p3cboxmodel.tex b/doc/handbook/ModelDescriptions/3p3cimplicitmodel.tex
similarity index 79%
rename from doc/handbook/ModelDescriptions/3p3cboxmodel.tex
rename to doc/handbook/ModelDescriptions/3p3cimplicitmodel.tex
index 91944a2ff6bd77243479ac517f63b719c2d46166..c8fddd9906af188aa5b40b5291bcf06839c1923b 100644
--- a/doc/handbook/ModelDescriptions/3p3cboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/3p3cimplicitmodel.tex
@@ -4,9 +4,9 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-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(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \]
+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} \mathbf{K} \left(\textbf{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_\alpha X_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_\alpha x_\alpha^\kappa \mbox{\bf K} (\textbf{grad}\, p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_\text{pm}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \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_\alpha X_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_\alpha x_\alpha^\kappa \mathbf{K} (\textbf{grad}\, p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_\text{pm}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\} - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha \end{eqnarray*}
 
 Note that these balance equations are molar.
 
@@ -14,7 +14,7 @@ All equations are discretized using a vertex-\/centered finite volume (box) or c
 
 The model uses commonly applied auxiliary conditions like $S_w + S_n + S_g = 1$ for the saturations and $x^w_\alpha + x^a_\alpha + x^c_\alpha = 1$ for the mole fractions. Furthermore, the phase pressures are related to each other via capillary pressures between the fluid phases, which are functions of the saturation, e.\-g. according to the approach of Parker et al.
 
-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\-:
+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.
diff --git a/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex b/doc/handbook/ModelDescriptions/3p3cniimplicitmodel.tex
similarity index 82%
rename from doc/handbook/ModelDescriptions/3p3cniboxmodel.tex
rename to doc/handbook/ModelDescriptions/3p3cniimplicitmodel.tex
index 87b28779f4df36296b60310fe1a9b2afb703308e..59743fb987cd93b9d9268609ff961e7bfe2e8fd5 100644
--- a/doc/handbook/ModelDescriptions/3p3cniboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/3p3cniimplicitmodel.tex
@@ -4,9 +4,9 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements three-\/phase three-\/component flow of three fluid phases $\alpha \in \{ \text{water, gas, NAPL} \}$ each composed of up to three components $\kappa \in \{ \text{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(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \]
+This model implements three-\/phase three-\/component flow of three fluid phases $\alpha \in \{ \text{water, gas, NAPL} \}$ each composed of up to three components $\kappa \in \{ \text{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} \mathbf{K} \left(\textbf{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_\alpha X_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_\alpha X_\alpha^\kappa \mbox{\bf K} (\textbf{grad}\; p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \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_\alpha X_\alpha^\kappa S_\alpha )}{\partial t} - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha} \varrho_\alpha X_\alpha^\kappa \mathbf{K} (\textbf{grad}\; p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\} \nonumber \\ \nonumber \\ && - \sum\limits_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\} - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha \end{eqnarray*}
 
 Note that these balance equations above are molar. In addition to that, a single balance of thermal energy is formulated for the fluid-\/filled porous medium under the assumption of local thermal equilibrium \begin{eqnarray*} && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\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 \mathbf{g} \right) \right\} \\ &-& \text{div} \left( \lambda_{pm} \textbf{grad} \, T \right) - q^h = 0 \qquad \alpha \in \{w, n, g\} \end{eqnarray*}
 
@@ -14,7 +14,7 @@ All equations are discretized using a vertex-\/centered finite volume (box) or c
 
 The model uses commonly applied auxiliary conditions like $S_w + S_n + S_g = 1$ for the saturations and $x^w_\alpha + x^a_\alpha + x^c_\alpha = 1$ for the mole fractions. Furthermore, the phase pressures are related to each other via capillary pressures between the fluid phases, which are functions of the saturation, e.\-g. according to the approach of Parker et al.
 
-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\-:
+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.
diff --git a/doc/handbook/ModelDescriptions/3pimplicitmodel.tex b/doc/handbook/ModelDescriptions/3pimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..1b78ef904bcaabc250657eb6cfff4ef46088c9b4
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/3pimplicitmodel.tex
@@ -0,0 +1,16 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements three-\/phase flow of three fluid phases $\alpha \in \{ water, gas, NAPL \}$ The standard multiphase Darcy approach is used as the equation for the conservation of momentum.
+
+By inserting this into the equations for the conservation of the components, the well-\/known multiphase flow equation is obtained.
+
+All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization.
+
+The model uses commonly applied auxiliary conditions like $S_w + S_n + S_g = 1$ for the saturations. Furthermore, the phase pressures are related to each other via capillary pressures between the fluid phases, which are functions of the saturation, e.\-g. according to the approach of Parker et al.
+
+The used primary variables are gas phase pressure $p_g$, water saturation $S_w$ and N\-A\-P\-L saturation $S_n$.
+
diff --git a/doc/handbook/ModelDescriptions/co2implicitmodel.tex b/doc/handbook/ModelDescriptions/co2implicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..80a8d27b742e3cd7d727102ef5b0fdb1fd3875e5
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/co2implicitmodel.tex
@@ -0,0 +1,8 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+See \hyperlink{a00557}{Two\-P\-Two\-C\-Model} for reference to the equations used. The \hyperlink{a00072}{C\-O2} model is derived from the 2p2c model. In the \hyperlink{a00072}{C\-O2} model the phase switch criterion is different from the 2p2c model. The phase switch occurs when the equilibrium concentration of a component in a phase is exceeded, instead of the sum of the components in the virtual phase (the phase which is not present) being greater that unity as done in the 2p2c model. The \hyperlink{a00079}{C\-O2\-Volume\-Variables} do not use a constraint solver for calculating the mole fractions as is the case in the 2p2c model. Instead mole fractions are calculated in the Fluid\-System with a given temperature, pressurem and salinity. The model is able to use either mole or mass fractions. The property use\-Moles can be set to either true or false in the problem file. Make sure that the according units are used in the problem setup. use\-Moles is set to false by default.
+
diff --git a/doc/handbook/ModelDescriptions/co2niimplicitmodel.tex b/doc/handbook/ModelDescriptions/co2niimplicitmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3cf083b448fe14e287d9d93309e18a6f44dc4a6d
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/co2niimplicitmodel.tex
@@ -0,0 +1,8 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+See Two\-P\-Two\-C\-N\-I model for reference to the equations. The C\-O2\-N\-I model is derived from the \hyperlink{a00072}{C\-O2} model. In the \hyperlink{a00072}{C\-O2} model the phase switch criterion is different from the 2p2c model. The phase switch occurs when the equilibrium concentration of a component in a phase is exceeded instead of the sum of the components in the virtual phase (the phase which is not present) being greater that unity as done in the 2p2c model. The \hyperlink{a00079}{C\-O2\-Volume\-Variables} do not use a constraint solver for calculating the mole fractions as is the case in the 2p2c model. Instead mole fractions are calculated in the Fluid\-System with a given temperature, pressure and salinity.
+
diff --git a/doc/handbook/ModelDescriptions/el1p2cmodel.tex b/doc/handbook/ModelDescriptions/el1p2cmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..23bdbac47efe11a0cac2b609b643107e9dbbe119
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/el1p2cmodel.tex
@@ -0,0 +1,32 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements a one-\/phase flow of an incompressible fluid, that consists of two components. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
+
+As an equation for the conservation of momentum within the fluid phase Darcy's approach is used\-: \[ v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right) \]
+
+Gravity can be enabled or disabled via the property system. By inserting this into the volume balance of the solid-\/fluid mixture, one gets \[ \frac{\partial \text{div} \textbf{u}}{\partial t} - \text{div} \left\{ \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right)\right\} = q \;, \]
+
+The transport of the components $\kappa \in \{ w, a \}$ is described by the following equation\-: \[ \frac{ \partial \phi_{eff} X^\kappa}{\partial t} - \text{div} \left\lbrace X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p - \varrho_w {\textbf g} \right) + D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa - \phi_{eff} X^\kappa \frac{\partial \boldsymbol{u}}{\partial t} \right\rbrace = q. \]
+
+If the model encounters stability problems, a stabilization term can be switched on. The stabilization term is defined in Aguilar et al (2008)\-: \[ \beta \text{div} \textbf{grad} \frac{\partial p}{\partial t} \] with $\beta$\-: \[ \beta = h^2 / 4(\lambda + 2 \mu) \] where $h$ is the discretization length.
+
+The balance equations with the stabilization term are given below\-: \[ \frac{\partial \text{div} \textbf{u}}{\partial t} - \text{div} \left\{ \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right) + \varrho_w \beta \textbf{grad} \frac{\partial p}{\partial t} \right\} = q \;, \]
+
+The transport of the components $\kappa \in \{ w, a \}$ is described by the following equation\-:
+
+\[ \frac{ \partial \phi_{eff} X^\kappa}{\partial t} - \text{div} \left\lbrace X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p - \varrho_w {\textbf g} \right) + \varrho_w X^\kappa \beta \textbf{grad} \frac{\partial p}{\partial t} + D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa - \phi_{eff} X^\kappa \frac{\partial \boldsymbol{u}}{\partial t} \right\rbrace = q. \]
+
+The quasi-\/stationary momentum balance equation is\-: \[ \text{div}\left( \boldsymbol{\sigma'}- p \boldsymbol{I} \right) + \left( \phi_{eff} \varrho_w + (1 - \phi_{eff}) * \varrho_s \right) {\textbf g} = 0 \;, \] with the effective stress\-: \[ \boldsymbol{\sigma'} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \boldsymbol{I}. \]
+
+and the strain tensor $\boldsymbol{\epsilon}$ as a function of the solid displacement gradient $\textbf{grad} \boldsymbol{u}$\-: \[ \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \boldsymbol{u} + \textbf{grad}^T \boldsymbol{u}). \]
+
+Here, the rock mechanics sign convention is switch off which means compressive stresses are $<$ 0 and tensile stresses are $>$ 0. The rock mechanics sign convention can be switched on for the vtk output via the property system.
+
+The effective porosity is calculated as a function of the solid displacement\-: \[ \phi_{eff} = \frac{\phi_{init} + \text{div} \boldsymbol{u}}{1 + \text{div}} \] All equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization.
+
+The primary variables are the pressure $p$ and the mole or mass fraction of dissolved component $x$ and the solid displacement vector $\boldsymbol{u}$.
+
diff --git a/doc/handbook/ModelDescriptions/el2pmodel.tex b/doc/handbook/ModelDescriptions/el2pmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..d95cdf204c4c11b30b2e60945220dc9e624cc939
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/el2pmodel.tex
@@ -0,0 +1,22 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements a two-\/phase flow of compressible immiscible fluids $\alpha \in \{ w, n \}$. The deformation of the solid matrix is described with a quasi-\/stationary momentum balance equation. The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941). The total stress acting on a rock is partially supported by the rock matrix and partially supported by the pore fluid. The effective stress represents the share of the total stress which is supported by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
+
+As an equation for the conservation of momentum within the fluid phases the standard multiphase Darcy's approach is used\-: \[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right) \]
+
+Gravity can be enabled or disabled via the property system. By inserting this into the continuity equation, one gets \[ \frac{\partial \phi_{eff} \varrho_\alpha S_\alpha}{\partial t} - \text{div} \left\{ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}_\text{eff} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right) - \phi_{eff} \varrho_\alpha S_\alpha \frac{\partial \mathbf{u}}{\partial t} \right\} - q_\alpha = 0 \;, \]
+
+A quasi-\/stationary momentum balance equation is solved for the changes with respect to the initial conditions (Darcis 2012), note that this implementation assumes the soil mechanics sign convention (i.\-e. compressive stresses are negative)\-: \[ \text{div}\left( \boldsymbol{\Delta \sigma'}- \Delta p_{eff} \boldsymbol{I} \right) + \Delta \varrho_b {\textbf g} = 0 \;, \] with the effective stress\-: \[ \boldsymbol{\sigma'} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \mathbf{I}. \]
+
+and the strain tensor $\boldsymbol{\epsilon}$ as a function of the solid displacement gradient $\textbf{grad} \mathbf{u}$\-: \[ \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \mathbf{u} + \textbf{grad}^T \mathbf{u}). \]
+
+Here, the rock mechanics sign convention is switch off which means compressive stresses are $<$ 0 and tensile stresses are $>$ 0. The rock mechanics sign convention can be switched on for the vtk output via the property system.
+
+The effective porosity and the effective permeability are calculated as a function of the solid displacement\-: \[ \phi_{eff} = \frac{\phi_{init} + \text{div} \mathbf{u}}{1 + \text{div} \mathbf{u}} \] \[ K_{eff} = K_{init} \text{exp}\left( 22.2(\phi_{eff}/\phi_{init} -1 )\right) \] The mass balance equations are discretized using a vertex-\/centered finite volume (box) or cell-\/centered finite volume scheme as spatial and the implicit Euler method as time discretization. The momentum balance equations are discretized using a standard Galerkin Finite Element method as spatial discretization scheme.
+
+The primary variables are the wetting phase pressure $p_w$, the nonwetting phase saturation $S_n$ and the solid displacement vector $\mathbf{u}$ (changes in solid displacement with respect to initial conditions).
+
diff --git a/doc/handbook/ModelDescriptions/elasticmodel.tex b/doc/handbook/ModelDescriptions/elasticmodel.tex
new file mode 100644
index 0000000000000000000000000000000000000000..79fa14958672f889b8afcb676cefbd59caf2fab5
--- /dev/null
+++ b/doc/handbook/ModelDescriptions/elasticmodel.tex
@@ -0,0 +1,14 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This file has been autogenerated from the LaTeX part of the   %
+% doxygen documentation; DO NOT EDIT IT! Change the model's .hh %
+% file instead!!                                                %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This model implements a linear elastic solid using Hooke's law as stress-\/strain relation and a quasi-\/stationary momentum balance equation\-: \[ \boldsymbol{\sigma} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \boldsymbol{I}. \]
+
+with the strain tensor $\boldsymbol{\epsilon}$ as a function of the solid displacement gradient $\textbf{grad} \boldsymbol{u}$\-: \[ \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \boldsymbol{u} + \textbf{grad}^T \boldsymbol{u}). \]
+
+Gravity can be enabled or disabled via the property system. By inserting this into the momentum balance equation, one gets \[ \text{div} \boldsymbol{\sigma} + \varrho {\textbf g} = 0 \;, \]
+
+The equation is discretized using a vertex-\/centered finite volume (box) scheme as spatial discretization.
+
diff --git a/doc/handbook/ModelDescriptions/mpncboxmodel.tex b/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
similarity index 96%
rename from doc/handbook/ModelDescriptions/mpncboxmodel.tex
rename to doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
index 35b633dde27293e22fe23bf795bc260e838763b1..bcdc19339173ea125b6e2e1ef2391582686b6a2b 100644
--- a/doc/handbook/ModelDescriptions/mpncboxmodel.tex
+++ b/doc/handbook/ModelDescriptions/mpncimplicitmodel.tex
@@ -6,7 +6,7 @@
 
 This model implements a $M$-\/phase flow of a fluid mixture composed of $N$ chemical species. The phases are denoted by lower index $\alpha \in \{ 1, \dots, M \}$. All fluid phases are mixtures of $N \geq M - 1$ chemical species which are denoted by the upper index $\kappa \in \{ 1, \dots, N \} $.
 
-The momentum approximation can be selected via \char`\"{}\-Base\-Flux\-Variables\char`\"{}\-: Darcy (\hyperlink{a00207}{Implicit\-Darcy\-Flux\-Variables}) and Forchheimer (\hyperlink{a00208}{Implicit\-Forchheimer\-Flux\-Variables}) relations are available for all Box models.
+The momentum approximation can be selected via \char`\"{}\-Base\-Flux\-Variables\char`\"{}\-: Darcy (\hyperlink{a00267}{Implicit\-Darcy\-Flux\-Variables}) and Forchheimer (\hyperlink{a00268}{Implicit\-Forchheimer\-Flux\-Variables}) relations are available for all Box models.
 
 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 \left(\varrho_\alpha x_\alpha^\kappa S_\alpha\right)}{\partial t} + \mathrm{div}\; \left\{ v_\alpha \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa \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 \]
 
diff --git a/doc/handbook/ModelDescriptions/richardsboxmodel.tex b/doc/handbook/ModelDescriptions/richardsimplicitmodel.tex
similarity index 100%
rename from doc/handbook/ModelDescriptions/richardsboxmodel.tex
rename to doc/handbook/ModelDescriptions/richardsimplicitmodel.tex
diff --git a/doc/handbook/ModelDescriptions/stokes2cmodel.tex b/doc/handbook/ModelDescriptions/stokes2cimplicitmodel.tex
similarity index 95%
rename from doc/handbook/ModelDescriptions/stokes2cmodel.tex
rename to doc/handbook/ModelDescriptions/stokes2cimplicitmodel.tex
index 9e63310fd12b0f59929f99b120c7b3fd96997f31..5ffdf2893a5c783a21eb134c03c04d598af1740a 100644
--- a/doc/handbook/ModelDescriptions/stokes2cmodel.tex
+++ b/doc/handbook/ModelDescriptions/stokes2cimplicitmodel.tex
@@ -10,7 +10,7 @@ Momentum Balance\-: \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}
 
 Mass balance equation\-: \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0 \]
 
-\hyperlink{a00070}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
+\hyperlink{a00085}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
 
 This is discretized using a fully-\/coupled vertex centered finite volume (box) scheme as spatial and the implicit Euler method as temporal discretization.
 
diff --git a/doc/handbook/ModelDescriptions/stokes2cnimodel.tex b/doc/handbook/ModelDescriptions/stokes2cniimplicitmodel.tex
similarity index 96%
rename from doc/handbook/ModelDescriptions/stokes2cnimodel.tex
rename to doc/handbook/ModelDescriptions/stokes2cniimplicitmodel.tex
index e84481ddf886b5127cac0a0f02ffc29ee6ed127f..7964642684c9109e09900e0c50004a7164259f7c 100644
--- a/doc/handbook/ModelDescriptions/stokes2cnimodel.tex
+++ b/doc/handbook/ModelDescriptions/stokes2cniimplicitmodel.tex
@@ -10,7 +10,7 @@ Momentum Balance\-: \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}
 
 Mass balance equation\-: \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0 \]
 
-\hyperlink{a00070}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
+\hyperlink{a00085}{Component} mass balance equation\-: \[ \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \boldsymbol{\nabla} x_g^\kappa \right) - q_g^\kappa = 0 \]
 
 Energy balance equation\-: \[ \frac{\partial (\varrho_g u_g)}{\partial t} + \boldsymbol{\nabla} \left( \boldsymbol{\cdot} \varrho_g h_g {\boldsymbol{v}}_g - \lambda_g \boldsymbol{\nabla} T \right) - q_T = 0 \]
 
diff --git a/doc/handbook/ModelDescriptions/stokesmodel.tex b/doc/handbook/ModelDescriptions/stokesimplicitmodel.tex
similarity index 50%
rename from doc/handbook/ModelDescriptions/stokesmodel.tex
rename to doc/handbook/ModelDescriptions/stokesimplicitmodel.tex
index faaf1fa4e9ca6dd05b5c060c83bc608019fc8469..fd999923e6f7f14be1e98c3f783f85bbb7e10940 100644
--- a/doc/handbook/ModelDescriptions/stokesmodel.tex
+++ b/doc/handbook/ModelDescriptions/stokesimplicitmodel.tex
@@ -4,11 +4,11 @@
 % file instead!!                                                %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-This model implements laminar Stokes flow of a single fluid, solving a momentum balance\-: \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(p_g {\bf {I}} - \mu_g \left(\boldsymbol{\nabla} \boldsymbol{v}_g + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)\right) - \varrho_g {\bf g} = 0, \]
+This model implements laminar Stokes flow of a single fluid, solving the momentum balance equation \[ \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(p_g {\bf {I}} - \mu_g \left(\boldsymbol{\nabla} \boldsymbol{v}_g + \boldsymbol{\nabla} \boldsymbol{v}_g^T\right)\right) - \varrho_g {\bf g} = 0, \]
 
-and the mass balance equation\-: \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0 \]
+and the mass balance equation \[ \frac{\partial \varrho_g}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0. \]
 
-By setting the property {\ttfamily Enable\-Navier\-Stokes} to {\ttfamily true} the Navier-\/\-Stokes equation can be solved. In this case an additional term is added to the momentum balance\-: \[ \varrho_g \left(\boldsymbol{v}_g \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{v}_g \]
+By setting the property {\ttfamily Enable\-Navier\-Stokes} to {\ttfamily true} the Navier-\/\-Stokes equation can be solved. In this case an additional term \[ \varrho_g \left(\boldsymbol{v}_g \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{v}_g \] is added to the momentum balance equation.
 
 This is discretized by a fully-\/coupled vertex-\/centered finite volume (box) scheme in space and by the implicit Euler method in time.
 
diff --git a/doc/handbook/models.tex b/doc/handbook/models.tex
index 6ca96fe34d8440c03db9e87078b459ebb8efae32..86b8d3b51d7601d10036a49b4d02bdd3575d9498 100644
--- a/doc/handbook/models.tex
+++ b/doc/handbook/models.tex
@@ -207,44 +207,65 @@ The fully-implicit models described in this section are using the box or the cel
 method as temporal discretization. The models themselves are located in
 subdirectories of \texttt{dumux/implicit} of the \Dumux distribution.
 
-\subsubsection{The single-phase model: OnePBoxModel} 
-\input{ModelDescriptions/1pboxmodel}
+\subsubsection{The single-phase model: OnePModel} 
+\input{ModelDescriptions/1pimplicitmodel}
 
-\subsubsection{The single-phase, two-component model:  OnePTwoCBoxModel} 
-\input{ModelDescriptions/1p2cboxmodel}
+\subsubsection{The single-phase, two-component model:  OnePTwoCModel} 
+\input{ModelDescriptions/1p2cimplicitmodel}
 
-\subsubsection{The two-phase model using the Richards assumption: RichardsBoxModel} 
-\input{ModelDescriptions/richardsboxmodel}
+\subsubsection{The two-phase model using the Richards assumption: RichardsModel} 
+\input{ModelDescriptions/richardsimplicitmodel}
 
-\subsubsection{The two-phase model: TwoPBoxModel}
-\input{ModelDescriptions/2pboxmodel}
+\subsubsection{The two-phase model: TwoPModel}
+\input{ModelDescriptions/2pimplicitmodel}
 
-\subsubsection{The non-isothermal two-phase model: TwoPNIBoxModel} 
-\input{ModelDescriptions/2pniboxmodel}
+\subsubsection{The non-isothermal two-phase model: TwoPNIModel} 
+\input{ModelDescriptions/2pniimplicitmodel}
 
-\subsubsection{The two-phase, two-component model: TwoPTwoCBoxModel} 
-\input{ModelDescriptions/2p2cboxmodel}
+\subsubsection{The two-phase, two-component model: TwoPTwoCModel} 
+\input{ModelDescriptions/2p2cimplicitmodel}
 
-\subsubsection{The non-isothermal two-phase, two-component model: TwoPTwoCNIBoxModel} 
-\input{ModelDescriptions/2p2cniboxmodel}
+\subsubsection{The CO2 model: CO2Model} 
+\input{ModelDescriptions/co2implicitmodel}
 
-\subsubsection{The three-phase, three-component model: ThreePThreeCBoxModel}
-\input{ModelDescriptions/3p3cboxmodel}
+\subsubsection{The non-isothermal two-phase, two-component model: TwoPTwoCNIModel} 
+\input{ModelDescriptions/2p2cniimplicitmodel}
 
-\subsubsection{The non-isothermal three-phase, three-component model: ThreePThreeCNIBoxModel} 
-\input{ModelDescriptions/3p3cniboxmodel}
+\subsubsection{The non-isothermal CO2 model: CO2NIModel} 
+\input{ModelDescriptions/co2niimplicitmodel}
 
-\subsubsection{The $M$-phase, $N$-component model: MpNcBoxModel} 
-\input{ModelDescriptions/mpncboxmodel}
+\subsubsection{The three-phase model: ThreePModel} 
+\input{ModelDescriptions/3pimplicitmodel}
+
+\subsubsection{The three-phase, three-component model: ThreePThreeCModel}
+\input{ModelDescriptions/3p3cimplicitmodel}
+
+\subsubsection{The non-isothermal three-phase, three-component model: ThreePThreeCNIModel} 
+\input{ModelDescriptions/3p3cniimplicitmodel}
+
+\subsubsection{The $M$-phase, $N$-component model: MpNcModel} 
+\input{ModelDescriptions/mpncimplicitmodel}
+
+\subsubsection{The two-phase, discrete fracture model: TwoPDFMModel} 
+\input{ModelDescriptions/2pdfmimplicitmodel}
 
 \subsubsection{The Stokes model: StokesModel} 
-\input{ModelDescriptions/stokesmodel}
+\input{ModelDescriptions/stokesimplicitmodel}
 
 \subsubsection{The isothermal two-component Stokes model: Stokes2cModel} 
-\input{ModelDescriptions/stokes2cmodel}
+\input{ModelDescriptions/stokes2cimplicitmodel}
 
 \subsubsection{The non-isothermal two-component Stokes model: Stokes2cniModel} 
-\input{ModelDescriptions/stokes2cnimodel}
+\input{ModelDescriptions/stokes2cniimplicitmodel}
+
+\subsubsection{The linear elastic model: ElasticModel} 
+\input{ModelDescriptions/elasticmodel}
+
+\subsubsection{The linear elastic one-phase two-component model: ElOnePTwoCModel} 
+\input{ModelDescriptions/el1p2cmodel}
+
+\subsubsection{The linear elastic two-phase model: ElTwoPModel} 
+\input{ModelDescriptions/el2pmodel}
 
 \subsection{Decoupled models}
 %
diff --git a/dumux/geomechanics/el2p/el2pmodel.hh b/dumux/geomechanics/el2p/el2pmodel.hh
index f3833e83e142bdab56ff7c0d12990bbc96842a98..8d6af6b4e0fe0f2ef9029fb6403e8830be889f4f 100644
--- a/dumux/geomechanics/el2p/el2pmodel.hh
+++ b/dumux/geomechanics/el2p/el2pmodel.hh
@@ -31,6 +31,12 @@
 #include <dune/pdelab/gridfunctionspace/interpolate.hh>
 
 namespace Dumux {
+
+namespace Properties {
+NEW_PROP_TAG(InitialDisplacement); //!< The initial displacement function
+NEW_PROP_TAG(InitialPressSat); //!< The initial pressure and saturation function
+}
+
 /*!
  * \ingroup ElTwoPBoxModel
  * \brief Adaption of the fully implicit scheme to the two-phase linear elasticity model.
@@ -93,12 +99,6 @@ namespace Dumux {
  * The primary variables are the wetting phase pressure \f$p_w\f$, the nonwetting phase saturation \f$S_n\f$ and the solid
  * displacement vector \f$\mathbf{u}\f$ (changes in solid displacement with respect to initial conditions).
  */
-
-namespace Properties {
-NEW_PROP_TAG(InitialDisplacement); //!< The initial displacement function
-NEW_PROP_TAG(InitialPressSat); //!< The initial pressure and saturation function
-}
-
 template<class TypeTag>
 class ElTwoPModel: public GET_PROP_TYPE(TypeTag, BaseModel)
 {