diff --git a/slides/img/model_domain_discretization.png b/slides/img/model_domain_discretization.png
new file mode 100644
index 0000000000000000000000000000000000000000..26db16c68e97ccb2c6f33c62e71a984fb2ed0c8e
Binary files /dev/null and b/slides/img/model_domain_discretization.png differ
diff --git a/slides/model.md b/slides/model.md
new file mode 100644
index 0000000000000000000000000000000000000000..b0157c4344556c07a61f3885bdc7e00c24259d43
--- /dev/null
+++ b/slides/model.md
@@ -0,0 +1,255 @@
+---
+title: Models in DuMux
+---
+
+# Implementing a new DuMu<sup>X</sup> model
+
+## What is a DuMu<sup>X</sup> model
+A DuMu<sup>X</sup> model is an implementation of a **mathematical model**, generally given by partial differential equations, by using **discretization** schemes.
+
+## What is a DuMu<sup>X</sup> model
+Mathematical model (PDE):
+$\begin{equation*}
+\small
+\frac{\partial S(u)}{\partial t} + \nabla \cdot \mathbf{F}(u) = q, \quad \forall (t,\mathbf{x}) \in (0,T] \times \Omega
+\end{equation*}$
+
+<img src=img/model_domain_discretization.png width="100%">
+
+Discrete model, e.g. using finite volumes:
+$\begin{equation*}
+\small
+|B| \frac{S_h(\mathbf{u}^{n+1}_h) - S_h(\mathbf{u}^{n}_h)}{\Delta t} + \sum_{\sigma \in \Sigma_B} F_{B,\sigma}(\mathbf{u}^{n+1}_h) = \int_{B} q^{n+1} \, dx, \quad \forall t_{n+1}\leq T, \; \forall B
+\end{equation*}$
+
+## What is a DuMu<sup>X</sup> model
+Discrete model, e.g. using finite volumes:
+$\begin{equation*}
+\small
+|B| \frac{S_h(\mathbf{u}^{n+1}_h) - S_h(\mathbf{u}^{n}_h)}{\Delta t} + \sum_{\sigma \in \Sigma_B} F_{B,\sigma}(\mathbf{u}^{n+1}_h) = \int_{B} q \, dx, \quad \forall t_{n+1}\leq T, \; \forall B
+\end{equation*}$
+
+* $S_h$: storage calculation
+* $F_{B,\sigma}$: flux over sub control volume face
+* $q$ source term
+
+How to implement these terms in DuMu<sup>X</sup>?
+
+**local residual**
+
+# Local Residual
+
+## Local Residual
+
+Implements these terms within
+
+* `computeStorage(...)`
+* `computeFlux(...)`
+* `computeSource(...)`
+
+
+## Example: Diffusion equation
+Mathematical model (PDE):
+
+$\begin{equation}
+\frac{\partial c}{\partial t} - \nabla \cdot (D \nabla c) = 0 \:\: \mathrm{in}\; \Omega \times (0,T] \\
+\end{equation}$
+
+with
+
+$\begin{aligned}
+ c: &\:\text{concentration} \\
+ D: &\:\text{constant diffusion coefficient} \\
+ \Omega: &\:\text{spatial domain} \\
+ T: &\:\text{end time}
+\end{aligned}$
+
+## Example: Diffusion equation
+Discrete model using the Box discretization:
+
+$\begin{equation}
+\vert B \vert \frac{c_B^{n+1}-c_B^{n}}{\Delta t} - \sum_{\sigma \in \Sigma_{B}} \left[ D \nabla c_h^{n+1} \cdot \boldsymbol{n}_{B,\sigma} \vert \sigma \vert \right] = 0,
+\end{equation}$
+
+with
+
+$\begin{aligned}
+ c_B^n: &\:\text{concentration at time $t_n$ and control volume $B$} \\
+ c_h: &\:\text{global discrete solution, interpolation using \textbf{basis functions}} \\
+ \mathbf{n}: &\:\text{unit outer normal vector} \\
+ \sigma: &\:\text{sub control volume face} \\
+\end{aligned}$
+
+
+## Local residual
+The local residual of the diffusion model:
+```cpp
+namespace Dumux {
+template<class TypeTag>
+class DiffusionModelLocalResidual
+: public GetPropType<TypeTag, Properties::BaseLocalResidual>
+{
+```
+
+## Local residual
+```cpp
+namespace Dumux {
+template<class TypeTag>
+class DiffusionModelLocalResidual
+: public GetPropType<TypeTag, Properties::BaseLocalResidual>
+{
+```
+Our class inherits from the class `BaseLocalResidual`, which has to be chosen depending on the discretization scheme, here *Box scheme*.
+
+## Storage term
+```cpp
+NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ NumEqVector storage;
+ storage[Indices::massBalanceEqIdx]
+ = volVars.priVar(Indices::concentrationIdx);
+ return storage;
+ }
+```
+
+## Flux term
+
+$\begin{equation}
+F_{B,\sigma} = -D \nabla c_h^{n+1} \cdot \boldsymbol{n}_{B,\sigma} \vert \sigma \vert
+\end{equation}$
+
+with
+
+$\begin{aligned}
+ c_h: &\:\text{discrete solution, interpolated using basis functions} \\
+ \sigma: &\:\text{sub control volume face} \\
+\end{aligned}$
+
+## Flux term
+Implementation takes place in the dedicated function `computeFlux`:
+```cpp
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ ...
+```
+
+## Flux term
+Implementation takes place in the dedicated function `computeFlux`:
+```cpp
+ NumEqVector computeFlux(...) const
+ {
+ // Compute ∇c
+ const auto& fluxVarCache = elemFluxVarsCache[scvf];
+ Dune::FieldVector<Scalar, dimWorld> gradConcentration(0.0);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ gradConcentration.axpy(
+ volVars.priVar(Indices::concentrationIdx),
+ fluxVarCache.gradN(scv.indexInElement())
+ );
+ }
+ ...
+ }
+```
+
+## Flux term
+Implementation takes place in the dedicated function `computeFlux`:
+```cpp
+ NumEqVector computeFlux(...) const
+ {
+ ...
+ NumEqVector flux;
+
+ // Compute the flux
+ flux[Indices::massBalanceEqIdx] = -1.0*vtmv(
+ scvf.unitOuterNormal(), problem.diffusionCoefficient(), gradConcentration
+ )*scvf.area();
+
+ return flux;
+
+ }
+```
+
+## Local Residual
+
+A *local residual* implements the discretized mathematical model.
+
+For its implementation different model-specific properties have to be set
+
+# Model properties
+
+## Model type tag
+
+The property tag is an empty struct with the respective name, e.g. `DiffusionModel`
+
+```cpp
+namespace Dumux::Properties::TTag {
+//! The diffusion model tag that we can specialize properties for
+struct DiffusionModel {};
+} // end namespace Dumux::Properties::TTag
+```
+## Model properties
+By specializing properties for the type tag `DiffusionModel`, every other class that knows about the type tag (this will be for example the assembler or the problem), can extract the type information that we specify here.
+
+All properties are defined within the namespace `Dumux::Properties`
+```cpp
+namespace Dumux::Properties {
+ //define all properties
+} // end namespace Dumux::Properties
+```
+
+## Defining model properties
+The type of the local residual is the class `DiffusionModelLocalResidual` defined from earlier
+```cpp
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::DiffusionModel>
+{ using type = DiffusionModelLocalResidual<TypeTag>; };
+```
+
+## Defining model properties
+The model traits specify information about the model:
+```cpp
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::DiffusionModel>
+{
+ struct type
+ {
+ struct Indices
+ {
+ static constexpr int concentrationIdx = 0;
+ static constexpr int massBalanceEqIdx = 0;
+ };
+ static constexpr int numEq() { return 1; }
+ };
+};
+```
+## Defining model properties
+Further model specific properties can be set accordingly by using the model property tag,
+i.e. `TTag::DiffusionModel`
+
+# Exercise: Model
+
+## Exercise: Model
+Implementation of a **nonlinear diffusion model**
+
+<figure>
+ <center>
+ <img src="../exercises/extradoc/exercisemodel_mri_denoise.gif" alt="denoising"/>
+ <figcaption> Denosing of MRI image using nonlinear diffusion model.</figcaption>
+ </center>
+</figure>
+
+## Tasks
+
+- Implement local residual
+- Set model properties
+- Use model in test case
+- Customize volume variables
\ No newline at end of file