diff --git a/examples/1ptracer/README.md b/examples/1ptracer/README.md
index d627b2a39567868531a094f0b440535039a7a863..21797824e612d6fa33564978ff697fdd45ba661f 100644
--- a/examples/1ptracer/README.md
+++ b/examples/1ptracer/README.md
@@ -40,13 +40,13 @@ are prescribed on a small strip close to the bottom boundary.
 ## Model description
 
 As mentioned above, two models are solved sequentially in this example. A single-phase
-model (_1p model_) is used to solve for the stationary velocity distribution of a fluid phase
-in the domain. The tracer transport is solved with the _tracer model_, which solves an advection-diffusion
+model _1p model_ ( [model documentation](https://dumux.org/docs/doxygen/master/group___one_p_model.html) ) is used to solve for the stationary velocity distribution of a fluid phase
+in the domain. The tracer transport is solved with the _tracer model_ ( [model documentation](https://dumux.org/docs/doxygen/master/group___tracer_model.html) ), which solves an advection-diffusion
 equation for a tracer component, which is assumed not to affect the density and viscosity
 of the fluid phase.
 
 ### 1p Model
-The single phase model uses Darcy's law as the equation for the momentum conservation:
+The single-phase flow model (model type tag `OneP`) uses Darcy's law as the equation for the momentum conservation:
 
 ```math
 \textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right),
@@ -70,7 +70,7 @@ in which both advective and diffusive transport mechanisms are considered:
 \phi \frac{ \partial \varrho X^\kappa}{\partial t} - \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0.
 ```
 
-Here, $`\textbf v`$ is a velocity field, which in this example is computed using the _1p model_ (see above). Moreover, $`X^\kappa`$ is the tracer mass fraction and $` D^\kappa_\text{pm} `$ is the
+Here, $`\textbf v`$ is the velocity field computed by the _1p model_ (see above). Moreover, $`X^\kappa`$ is the tracer mass fraction and $` D^\kappa_\text{pm} `$ is the
 effective diffusivity. In this example, the effective diffusivity is a function of the diffusion
 coefficient of the tracer component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous
 medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
@@ -84,8 +84,8 @@ The primary variable used in this model is the tracer mass fraction $`X^\kappa`$
 ### Discretization
 
 In this example, all equations are discretized using cell-centered finite volumes with two-point flux
-approximation as spatial discretization scheme. For details on the discretization schemes available in
-DuMuX, have a look at the [code documentation](https://dumux.org/docs/doxygen/master/group___discretization.html).
+approximation as spatial discretization scheme (`CCTpfaModel`). For details on the discretization schemes available in
+DuMu<sup>x</sup>, have a look at the [code documentation](https://dumux.org/docs/doxygen/master/group___discretization.html).
 We use the [implicit Euler method](https://dumux.org/docs/doxygen/master/basic-numerics.html) as
 time discretization scheme for the tracer component balance equation solved in the _tracer model_.
 
@@ -98,7 +98,7 @@ In the following, we take a closer look at the source files for this example:
     ├── CMakeLists.txt          -> build system file
     ├── main.cc                 -> main program flow
     ├── params.input            -> runtime parameters
-    ├── properties._1p.hh       -> compile time settings for the single-phase flow simulation
+    ├── properties_1p.hh       -> compile time settings for the single-phase flow simulation
     ├── problem_1p.hh           -> boundary & initial conditions for the single-phase flow simulation
     ├── spatialparams_1p.hh     -> parameter distributions for the single-phase flow simulation
     ├── properties_tracer.hh    -> compile time settings for the tracer transport simulation
@@ -106,8 +106,8 @@ In the following, we take a closer look at the source files for this example:
     └── spatialparams_tracer.hh -> parameter distributions for the tracer transport simulation
 ```
 
-In order to define a simulation setup in DuMuX, you need to implement compile-time settings,
-where you specify the classes and compile-time options that DuMuX should use for the simulation.
+In order to define a simulation setup in DuMu<sup>x</sup>, you need to implement compile-time settings,
+where you specify the classes and compile-time options that DuMu<sup>x</sup> should use for the simulation.
 Moreover, a `Problem` class needs to be implemented, in which the initial and boundary conditions
 are specified. Finally, spatially-distributed values for the parameters required by the used model
 are implemented in a `SpatialParams` class.
diff --git a/examples/1ptracer/doc/1p.md b/examples/1ptracer/doc/1p.md
index a996628c72551330985f0caa0aa2ec58eaf9f210..564a97dde9e613d1dd418008ba73085cd7536d93 100644
--- a/examples/1ptracer/doc/1p.md
+++ b/examples/1ptracer/doc/1p.md
@@ -9,7 +9,7 @@
 The single-phase flow setup is implemented in the files `properties_1p.hh`,
 `problem_1p.hh` and `spatialparams_1p.hh`. In the first of these files, a new
 type tag is declared for this problem. This then allows the specialization
-of DuMuX `properties` for this type tag, which can be used to customize compile-time
+of DuMu<sup>x</sup> `properties` for this type tag, which can be used to customize compile-time
 settings for the simulation. Two exemplary `properties`, that are mandatory to be
 specialized, are `Problem` and `SpatialParams`. With the first, one sets the
 `Problem` class to be used, in which users can define initial and boundary conditions.
@@ -35,7 +35,7 @@ for which we then specialize `properties` to the needs of the desired setup.
 <details><summary> Click to show includes</summary>
 
 The `OneP` type tag specializes most of the `properties` required for single-
-phase flow simulations in DuMuX. We will use this in the following to inherit the
+phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
 respective properties and subsequently specialize those `properties` for our
 type tag, which we want to modify or for which no meaningful default can be set.
 
@@ -137,7 +137,7 @@ public:
 
 
 These are all mandatory `properties`. However, we also want to specialize the
-`LocalResidual` property, as we want to use analytic differentation (see `main.cc`).
+`LocalResidual` property, as we want to use analytic differentiation (see `main.cc`).
 The default for the `LocalResidual` property, specialized by the type tag `OneP`,
 is the implementation of the local residual for compressible single-phase flow.
 Therein, the functionality required for analytic derivatives is not implemented
@@ -155,7 +155,7 @@ Apart from that, we also set some properties related to memory management
 throughout the simulation.
 <details><summary> Click to show caching properties</summary>
 
-In Dumux, one has the option to activate/deactivate the grid-wide caching of
+In DuMu<sup>x<\sup>, one has the option to activate/deactivate the grid-wide caching of
 geometries and variables. If active, the CPU time can be significantly reduced
 as less dynamic memory allocation procedures are necessary. Per default, grid-wide
 caching is disabled to ensure minimal memory requirements, however, in this example we
@@ -209,7 +209,7 @@ Include the `BoundaryTypes` class which specifies the boundary types set in this
 ```
 
 ### The problem class
-We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+We enter the problem class where all necessary boundary and initial conditions are set for our simulation.
 As we are solving a problem related to flow in porous media, we inherit from the base class `PorousMediumFlowProblem`.
 
 ```cpp
@@ -252,7 +252,7 @@ equations on the same boundary) are not accepted for cell-centered finite volume
         // we define a small epsilon value
         Scalar eps = 1.0e-6;
 
-        // Initially, set Neumann boundary conditions for all equations
+        // Initially, set Neumann boundary conditions for all equations.
         BoundaryTypes values;
         values.setAllNeumann();
 
@@ -272,7 +272,7 @@ set a pressure of 1.1 bar and 1 bar at the bottom and top boundaries, respective
 ```cpp
     PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
     {
-        // instantiate a primary variables object
+        // Instantiate a primary variables object
         PrimaryVariables values;
 
         // and assign a pressure value in [Pa] such that at the bottom boundary
@@ -280,16 +280,11 @@ set a pressure of 1.1 bar and 1 bar at the bottom and top boundaries, respective
         values[0] = 1.0e5*(1.1 - globalPos[dimWorld-1]*0.1);
         return values;
     }
-```
-
-
-```cpp
 
 }; // end class definition of OnePTestProblem
 } // end namespace Dumux
 ```
 
-[[/codeblock]]
 
 </details>
 
@@ -299,7 +294,7 @@ set a pressure of 1.1 bar and 1 bar at the bottom and top boundaries, respective
 
 This file contains the __spatial parameters class__ which defines the
 distributions for the porous medium parameters permeability and porosity
-over the computational grid
+over the computational grid.
 
 
 <details open>
@@ -310,7 +305,7 @@ over the computational grid
 In this example, we use a randomly generated and element-wise distributed
 permeability field. For this, we use the random number generation facilities
 provided by the C++ standard library.
-Dumux additionally implements some simpler but
+DuMu<sup>x</sup> additionally implements some simpler but
 standard-library-implementation-independent random distributions
 that are accurate enough for our purposes
 but allow to generate reproducible and portable random fields (when using a constant seed)
@@ -333,7 +328,7 @@ example will inherit.
 ### The spatial parameters class
 
 In the `OnePTestSpatialParams` class, we define all functions needed to describe
-the porous medium, e.g. porosity and permeability, for the single-phase problem.
+the porous medium, e.g. porosity and permeability, for the simulated single-phase flow scenario.
 We inherit from the `FVPorousMediumFlowSpatialParamsOneP` class here, which is the base class
 for spatial parameters in the context of single-phase porous medium flow
 applications using finite volume discretization schemes.
@@ -346,7 +341,7 @@ class OnePTestSpatialParams
 : public FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar,
                              OnePTestSpatialParams<GridGeometry, Scalar>>
 {
-    // The following convenience aliases will be used throughout this class
+    // The following convenience aliases will be used throughout this class.
     using GridView = typename GridGeometry::GridView;
     using FVElementGeometry = typename GridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -369,7 +364,7 @@ We generate a random permeability field upon construction of the spatial paramet
 using lognormal distributions. The mean values for the permeability inside and outside of a
 low-permeable lens (given by the coordinates `lensLowerLeft_` and `lensUpperRight_`) are defined
 in the variables  `permeabilityLens_` and `permeability_`. The respective values are obtained
-from the input file making use of the free function `getParam`. We use a standard deviarion
+from the input file (`params.input`) making use of the free function `getParam`. We use a standard deviarion
 of 10% here and compute permeabily values for all elements of the computational grid.
 
 ```cpp
diff --git a/examples/1ptracer/doc/1p_intro.md b/examples/1ptracer/doc/1p_intro.md
index 8ba21fbc1e0168672a70253def48f5682b64647e..fe6f5cee9ac230d3f9ef4a4aaaf86bfed0b0eea6 100644
--- a/examples/1ptracer/doc/1p_intro.md
+++ b/examples/1ptracer/doc/1p_intro.md
@@ -3,7 +3,7 @@
 The single-phase flow setup is implemented in the files `properties_1p.hh`,
 `problem_1p.hh` and `spatialparams_1p.hh`. In the first of these files, a new
 type tag is declared for this problem. This then allows the specialization
-of DuMuX `properties` for this type tag, which can be used to customize compile-time
+of DuMu<sup>x</sup> `properties` for this type tag, which can be used to customize compile-time
 settings for the simulation. Two exemplary `properties`, that are mandatory to be
 specialized, are `Problem` and `SpatialParams`. With the first, one sets the
 `Problem` class to be used, in which users can define initial and boundary conditions.
diff --git a/examples/1ptracer/doc/_intro.md b/examples/1ptracer/doc/_intro.md
index 9aeb98ff7fa862f71b4702e4ec65e7e587780dcb..d761b601a943f0e9baa09a550456f1b7ef2db1d0 100644
--- a/examples/1ptracer/doc/_intro.md
+++ b/examples/1ptracer/doc/_intro.md
@@ -38,13 +38,13 @@ are prescribed on a small strip close to the bottom boundary.
 ## Model description
 
 As mentioned above, two models are solved sequentially in this example. A single-phase
-model (_1p model_) is used to solve for the stationary velocity distribution of a fluid phase
-in the domain. The tracer transport is solved with the _tracer model_, which solves an advection-diffusion
+model _1p model_ ( [model documentation](https://dumux.org/docs/doxygen/master/group___one_p_model.html) ) is used to solve for the stationary velocity distribution of a fluid phase
+in the domain. The tracer transport is solved with the _tracer model_ ( [model documentation](https://dumux.org/docs/doxygen/master/group___tracer_model.html) ), which solves an advection-diffusion
 equation for a tracer component, which is assumed not to affect the density and viscosity
 of the fluid phase.
 
 ### 1p Model
-The single phase model uses Darcy's law as the equation for the momentum conservation:
+The single-phase flow model (model type tag `OneP`) uses Darcy's law as the equation for the momentum conservation:
 
 ```math
 \textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right),
@@ -68,7 +68,7 @@ in which both advective and diffusive transport mechanisms are considered:
 \phi \frac{ \partial \varrho X^\kappa}{\partial t} - \text{div} \left\lbrace \varrho X^\kappa {\textbf v} + \varrho D^\kappa_\text{pm} \textbf{grad} X^\kappa \right\rbrace = 0.
 ```
 
-Here, $`\textbf v`$ is a velocity field, which in this example is computed using the _1p model_ (see above). Moreover, $`X^\kappa`$ is the tracer mass fraction and $` D^\kappa_\text{pm} `$ is the
+Here, $`\textbf v`$ is the velocity field computed by the _1p model_ (see above). Moreover, $`X^\kappa`$ is the tracer mass fraction and $` D^\kappa_\text{pm} `$ is the
 effective diffusivity. In this example, the effective diffusivity is a function of the diffusion
 coefficient of the tracer component $`D^\kappa`$ and the porosity and tortuosity $`\tau`$ of the porous
 medium (see [dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh)):
@@ -82,8 +82,8 @@ The primary variable used in this model is the tracer mass fraction $`X^\kappa`$
 ### Discretization
 
 In this example, all equations are discretized using cell-centered finite volumes with two-point flux
-approximation as spatial discretization scheme. For details on the discretization schemes available in
-DuMuX, have a look at the [code documentation](https://dumux.org/docs/doxygen/master/group___discretization.html).
+approximation as spatial discretization scheme (`CCTpfaModel`). For details on the discretization schemes available in
+DuMu<sup>x</sup>, have a look at the [code documentation](https://dumux.org/docs/doxygen/master/group___discretization.html).
 We use the [implicit Euler method](https://dumux.org/docs/doxygen/master/basic-numerics.html) as
 time discretization scheme for the tracer component balance equation solved in the _tracer model_.
 
@@ -96,7 +96,7 @@ In the following, we take a closer look at the source files for this example:
     ├── CMakeLists.txt          -> build system file
     ├── main.cc                 -> main program flow
     ├── params.input            -> runtime parameters
-    ├── properties._1p.hh       -> compile time settings for the single-phase flow simulation
+    ├── properties_1p.hh       -> compile time settings for the single-phase flow simulation
     ├── problem_1p.hh           -> boundary & initial conditions for the single-phase flow simulation
     ├── spatialparams_1p.hh     -> parameter distributions for the single-phase flow simulation
     ├── properties_tracer.hh    -> compile time settings for the tracer transport simulation
@@ -104,8 +104,8 @@ In the following, we take a closer look at the source files for this example:
     └── spatialparams_tracer.hh -> parameter distributions for the tracer transport simulation
 ```
 
-In order to define a simulation setup in DuMuX, you need to implement compile-time settings,
-where you specify the classes and compile-time options that DuMuX should use for the simulation.
+In order to define a simulation setup in DuMu<sup>x</sup>, you need to implement compile-time settings,
+where you specify the classes and compile-time options that DuMu<sup>x</sup> should use for the simulation.
 Moreover, a `Problem` class needs to be implemented, in which the initial and boundary conditions
 are specified. Finally, spatially-distributed values for the parameters required by the used model
 are implemented in a `SpatialParams` class.
diff --git a/examples/1ptracer/doc/main.md b/examples/1ptracer/doc/main.md
index a66dcf25fccc7c8d143f703ad111117c6511b6e4..74b1bc9b7eb081bd07572343bf85b6854ab38d5f 100644
--- a/examples/1ptracer/doc/main.md
+++ b/examples/1ptracer/doc/main.md
@@ -52,7 +52,7 @@ systems arising from finite volume discretizations (box-scheme, tpfa-approximati
 #include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
 ```
 
-The following class provides a convenient way of writing of dumux simulation results to VTK format.
+The following class provides a convenient way of writing DuMu<sup>x</sup> simulation results to VTK format.
 
 ```cpp
 #include <dumux/io/vtkoutputmodule.hh>
diff --git a/examples/1ptracer/doc/tracer.md b/examples/1ptracer/doc/tracer.md
index a69f7f573266e1f50aa0cefec666f07697697ee7..3159fa40f890d2576623d067f8cebab0d4490884 100644
--- a/examples/1ptracer/doc/tracer.md
+++ b/examples/1ptracer/doc/tracer.md
@@ -9,7 +9,7 @@
 The tracer transport setup is implemented in the files `properties_tracer.hh`,
 `problem_tracer.hh` and `spatialparams_tracer.hh`. In the first of these files, a new
 `TypeTag` is declared for this problem. This then allows the specialization
-of DuMuX `properties` for this new `TypeTag`, which can be used to customize compile-time
+of DuMu<sup>x</sup> `properties` for this new `TypeTag`, which can be used to customize compile-time
 settings for the simulation. Two exemplary `properties`, that are mandatory to be
 specialized, are `Problem` and `SpatialParams`. With the first, one sets the
 `Problem` class to be used, in which users can define initial and boundary conditions.
@@ -35,7 +35,7 @@ which we then specialize `properties` to the needs of the desired setup.
 <details><summary> Click to show includes</summary>
 As for the single-phase problem, atype tag is defined also for this simulation.
 Here, we inherit all properties of the `Tracer` type tag, a convenience type tag
-that specializes most of the required properties for tracer transport flow simulations in DuMuX.
+that specializes most of the required properties for tracer transport flow simulations in DuMu<sup>x</sup>.
 
 ```cpp
 #include <dumux/porousmediumflow/tracer/model.hh>
diff --git a/examples/1ptracer/doc/tracer_intro.md b/examples/1ptracer/doc/tracer_intro.md
index f901597d6691974e3cc2cef9c097dcf5dff963e2..e6d9f4d063911d3a8787b115bcf428f727750917 100644
--- a/examples/1ptracer/doc/tracer_intro.md
+++ b/examples/1ptracer/doc/tracer_intro.md
@@ -3,7 +3,7 @@
 The tracer transport setup is implemented in the files `properties_tracer.hh`,
 `problem_tracer.hh` and `spatialparams_tracer.hh`. In the first of these files, a new
 `TypeTag` is declared for this problem. This then allows the specialization
-of DuMuX `properties` for this new `TypeTag`, which can be used to customize compile-time
+of DuMu<sup>x</sup> `properties` for this new `TypeTag`, which can be used to customize compile-time
 settings for the simulation. Two exemplary `properties`, that are mandatory to be
 specialized, are `Problem` and `SpatialParams`. With the first, one sets the
 `Problem` class to be used, in which users can define initial and boundary conditions.
diff --git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc
index 1482d210491c70b3eb75d237ebdb924dbf5f6f24..8de03e6172f23611bccce1c0de31da629f441668 100644
--- a/examples/1ptracer/main.cc
+++ b/examples/1ptracer/main.cc
@@ -33,7 +33,7 @@
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
 
-// The following class provides a convenient way of writing of dumux simulation results to VTK format.
+// The following class provides a convenient way of writing DuMu<sup>x</sup> simulation results to VTK format.
 #include <dumux/io/vtkoutputmodule.hh>
 // The gridmanager constructs a grid from the information in the input or grid file.
 // Many different Dune grid implementations are supported, of which a list can be found
diff --git a/examples/1ptracer/problem_1p.hh b/examples/1ptracer/problem_1p.hh
index 2a0535f543571c35dc1ac9067628b2e87a5a97eb..eea7f96c9f99747146cc4941c33a9a1384059988 100644
--- a/examples/1ptracer/problem_1p.hh
+++ b/examples/1ptracer/problem_1p.hh
@@ -24,7 +24,7 @@
 #include <dumux/common/boundarytypes.hh>
 
 // ### The problem class
-// We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+// We enter the problem class where all necessary boundary and initial conditions are set for our simulation.
 // As we are solving a problem related to flow in porous media, we inherit from the base class `PorousMediumFlowProblem`.
 // [[codeblock]]
 namespace Dumux {
@@ -65,7 +65,7 @@ public:
         // we define a small epsilon value
         Scalar eps = 1.0e-6;
 
-        // Initially, set Neumann boundary conditions for all equations
+        // Initially, set Neumann boundary conditions for all equations.
         BoundaryTypes values;
         values.setAllNeumann();
 
@@ -84,7 +84,7 @@ public:
     // [[codeblock]]
     PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
     {
-        // instantiate a primary variables object
+        // Instantiate a primary variables object
         PrimaryVariables values;
 
         // and assign a pressure value in [Pa] such that at the bottom boundary
@@ -92,7 +92,6 @@ public:
         values[0] = 1.0e5*(1.1 - globalPos[dimWorld-1]*0.1);
         return values;
     }
-    // [[/codeblock]]
 
 }; // end class definition of OnePTestProblem
 } // end namespace Dumux
diff --git a/examples/1ptracer/properties_1p.hh b/examples/1ptracer/properties_1p.hh
index 51a943f59add504e8c776591e00ee6fa2b82fc48..ccb5d8dd55001672deaef9c2338d86779ea6c537 100644
--- a/examples/1ptracer/properties_1p.hh
+++ b/examples/1ptracer/properties_1p.hh
@@ -19,7 +19,7 @@
 // [[details]] includes
 //
 // The `OneP` type tag specializes most of the `properties` required for single-
-// phase flow simulations in DuMuX. We will use this in the following to inherit the
+// phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
 // respective properties and subsequently specialize those `properties` for our
 // type tag, which we want to modify or for which no meaningful default can be set.
 #include <dumux/porousmediumflow/1p/model.hh>
@@ -97,7 +97,7 @@ public:
 // [[/codeblock]]
 //
 // These are all mandatory `properties`. However, we also want to specialize the
-// `LocalResidual` property, as we want to use analytic differentation (see `main.cc`).
+// `LocalResidual` property, as we want to use analytic differentiation (see `main.cc`).
 // The default for the `LocalResidual` property, specialized by the type tag `OneP`,
 // is the implementation of the local residual for compressible single-phase flow.
 // Therein, the functionality required for analytic derivatives is not implemented
@@ -112,7 +112,7 @@ struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncom
 // throughout the simulation.
 // [[details]] caching properties
 //
-// In Dumux, one has the option to activate/deactivate the grid-wide caching of
+// In DuMu<sup>x<\sup>, one has the option to activate/deactivate the grid-wide caching of
 // geometries and variables. If active, the CPU time can be significantly reduced
 // as less dynamic memory allocation procedures are necessary. Per default, grid-wide
 // caching is disabled to ensure minimal memory requirements, however, in this example we
diff --git a/examples/1ptracer/properties_tracer.hh b/examples/1ptracer/properties_tracer.hh
index bfe863aec3188cc542197e98ac9da4b52ed9f979..dd091c04e3e28d0a6a0ad7eadb841ab0471c659f 100644
--- a/examples/1ptracer/properties_tracer.hh
+++ b/examples/1ptracer/properties_tracer.hh
@@ -19,7 +19,7 @@
 // [[details]] includes
 // As for the single-phase problem, atype tag is defined also for this simulation.
 // Here, we inherit all properties of the `Tracer` type tag, a convenience type tag
-// that specializes most of the required properties for tracer transport flow simulations in DuMuX.
+// that specializes most of the required properties for tracer transport flow simulations in DuMu<sup>x</sup>.
 #include <dumux/porousmediumflow/tracer/model.hh>
 
 // We use YaspGrid, an implementation of the dune grid interface for structured grids.
diff --git a/examples/1ptracer/spatialparams_1p.hh b/examples/1ptracer/spatialparams_1p.hh
index 40091ec3f8f184617f1b9e4de03885dd05729401..86a34ece117dd5c668e2471518a082b6616de8fd 100644
--- a/examples/1ptracer/spatialparams_1p.hh
+++ b/examples/1ptracer/spatialparams_1p.hh
@@ -12,7 +12,7 @@
 //
 // This file contains the __spatial parameters class__ which defines the
 // distributions for the porous medium parameters permeability and porosity
-// over the computational grid
+// over the computational grid.
 //
 // [[content]]
 //
@@ -20,7 +20,7 @@
 // In this example, we use a randomly generated and element-wise distributed
 // permeability field. For this, we use the random number generation facilities
 // provided by the C++ standard library.
-// Dumux additionally implements some simpler but
+// DuMu<sup>x</sup> additionally implements some simpler but
 // standard-library-implementation-independent random distributions
 // that are accurate enough for our purposes
 // but allow to generate reproducible and portable random fields (when using a constant seed)
@@ -36,7 +36,7 @@
 // ### The spatial parameters class
 //
 // In the `OnePTestSpatialParams` class, we define all functions needed to describe
-// the porous medium, e.g. porosity and permeability, for the single-phase problem.
+// the porous medium, e.g. porosity and permeability, for the simulated single-phase flow scenario.
 // We inherit from the `FVPorousMediumFlowSpatialParamsOneP` class here, which is the base class
 // for spatial parameters in the context of single-phase porous medium flow
 // applications using finite volume discretization schemes.
@@ -48,7 +48,7 @@ class OnePTestSpatialParams
 : public FVPorousMediumFlowSpatialParamsOneP<GridGeometry, Scalar,
                              OnePTestSpatialParams<GridGeometry, Scalar>>
 {
-    // The following convenience aliases will be used throughout this class
+    // The following convenience aliases will be used throughout this class.
     using GridView = typename GridGeometry::GridView;
     using FVElementGeometry = typename GridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -70,7 +70,7 @@ public:
     // using lognormal distributions. The mean values for the permeability inside and outside of a
     // low-permeable lens (given by the coordinates `lensLowerLeft_` and `lensUpperRight_`) are defined
     // in the variables  `permeabilityLens_` and `permeability_`. The respective values are obtained
-    // from the input file making use of the free function `getParam`. We use a standard deviarion
+    // from the input file (`params.input`) making use of the free function `getParam`. We use a standard deviarion
     // of 10% here and compute permeabily values for all elements of the computational grid.
     // [[codeblock]]
     OnePTestSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
diff --git a/examples/2pinfiltration/README.md b/examples/2pinfiltration/README.md
index 20da8e72f56b837e8e9f1fa88265b23e32d5f64a..002dff365b2e52d5f3892e697233dea103e65e01 100644
--- a/examples/2pinfiltration/README.md
+++ b/examples/2pinfiltration/README.md
@@ -45,7 +45,7 @@ Van Genuchten-Mualem relationships (see
 [Van Genuchten (1980)](https://doi.org/10.2136/sssaj1980.03615995004400050002x)
 and
 [Mualem (1976)](https://doi.org/10.1029/WR012i003p00513))
-for the capillary pressure $`pc = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
+for the capillary pressure $`p_c = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
 The parameters for these relationships are specified in the `spatialparams.hh` file.
 
 With the additional constraint that $`S_w + S_n = 1`$, the number of unknowns is reduced to two.
diff --git a/examples/2pinfiltration/doc/2p.md b/examples/2pinfiltration/doc/2p.md
index 5f522154db4a36ed0d315bd66b6a0e2292d192bc..111bb72030f28a618a735f1c5c4cc5eb1d82928c 100644
--- a/examples/2pinfiltration/doc/2p.md
+++ b/examples/2pinfiltration/doc/2p.md
@@ -93,7 +93,7 @@ struct SpatialParams<TypeTag, TTag::PointSourceExample>
 The `Grid` property tells the
 simulator to use ALUGrid - an unstructured grid manager - here
 configured for grid and coordinate dimensions `2`,
-hexahedral element types (`Dune::cube`) and non-conforming refinement mode.
+hexahedral element types (`Dune::cube`) and non-conforming refinement mode (`Dune::nonconforming`).
 `Dune::ALUGrid` is declared in the included header `dune/alugrid/grid.hh`
 from the Dune module `dune-alugrid`.
 
@@ -274,8 +274,8 @@ and the DNAPL saturation.
 
 #### Neumann boundaries
 In our case, we need to specify mass fluxes for our two liquid phases.
-Negative sign means influx and the unit of the boundary flux is $`kg/(m^2 s)`$.
-On the inlet area, we set a DNAPL influx of $`0.04 kg/(m^2 s)`$. On all other
+Negative sign means influx and the unit of the boundary flux is $\mathrm{kg}/\mathrm{(m^2 s)}$.
+On the inlet area, we set a DNAPL influx of $0.04 \, \mathrm{kg}/\mathrm{(m^2 s)}$. On all other
 Neumann boundaries, the boundary flux is zero.
 
 ```cpp
@@ -313,7 +313,7 @@ In this scenario, we set a point source (e.g. modeling a well). The point source
 Point sources are added by pushing them into the vector `pointSources`.
 The `PointSource` constructor takes two arguments.
 The first argument is a coordinate array containing the position in space,
-the second argument is an array of source value for each equation (in units of $`kg/s`$).
+the second argument is an array of source values for each equation (in units of $`\mathrm{kg}/\mathrm{s}`$).
 Recall that the first equation is the water phase mass balance
 and the second equation is the DNAPL phase mass balance.
 
diff --git a/examples/2pinfiltration/doc/_intro.md b/examples/2pinfiltration/doc/_intro.md
index 75982c3351463060846f5ebe096eb0a17cae11d1..a3705132df17f62e2a4a4cfa46064e423d4098ea 100644
--- a/examples/2pinfiltration/doc/_intro.md
+++ b/examples/2pinfiltration/doc/_intro.md
@@ -43,7 +43,7 @@ Van Genuchten-Mualem relationships (see
 [Van Genuchten (1980)](https://doi.org/10.2136/sssaj1980.03615995004400050002x)
 and
 [Mualem (1976)](https://doi.org/10.1029/WR012i003p00513))
-for the capillary pressure $`pc = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
+for the capillary pressure $`p_c = p_n - p_w`$ and the relative permeabilities $`k_r\alpha`$.
 The parameters for these relationships are specified in the `spatialparams.hh` file.
 
 With the additional constraint that $`S_w + S_n = 1`$, the number of unknowns is reduced to two.
diff --git a/examples/2pinfiltration/doc/main.md b/examples/2pinfiltration/doc/main.md
index 468bbb09c7ed72c0f3080e796fa07f222da21d03..e2a02bb8c18a6b5455be3d59d48e7018bcb9c74c 100644
--- a/examples/2pinfiltration/doc/main.md
+++ b/examples/2pinfiltration/doc/main.md
@@ -20,7 +20,7 @@ The code documentation is structured as follows:
 <summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
 
 
-This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
+This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method.
 ### Included header files
 <details><summary> Click to show includes</summary>
 
@@ -45,7 +45,7 @@ In Dumux, a property system is used to specify the model. For this, different pr
 #include <dumux/common/initialize.hh>
 ```
 
-We include the linear solver to be used to solve the linear system and the nonlinear  Newton's method
+We include the linear solver to be used to solve the linear system and the nonlinear Newton's method.
 
 ```cpp
 #include <dumux/linear/istlsolvers.hh>
@@ -54,7 +54,7 @@ We include the linear solver to be used to solve the linear system and the nonli
 #include <dumux/nonlinear/newtonsolver.hh>
 ```
 
-Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a file that defines the different differentiation methods used to compute the derivatives of the residual
+Further, we include the assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a header file that defines the different differentiation methods used to compute the derivatives of the residual.
 
 ```cpp
 #include <dumux/assembly/fvassembler.hh>
@@ -73,7 +73,7 @@ The gridmanager constructs a grid from the information in the input or grid file
 #include <dumux/io/grid/gridmanager_alu.hh>
 ```
 
-We include several files which are needed for the adaptive grid
+We include several files which are needed for the adaptive grid.
 
 ```cpp
 #include <dumux/adaptive/adapt.hh>
@@ -83,7 +83,7 @@ We include several files which are needed for the adaptive grid
 #include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
 ```
 
-Finally, we include the properties which configure the simulation
+Finally, we include the properties which configure the simulation.
 
 ```cpp
 #include "properties.hh"
@@ -122,13 +122,13 @@ file or specified domain dimensions and number of cells, as done in this example
 
     // The instationary non-linear problem is run on this grid.
     //
-    // we compute on the leaf grid view
+    // We compute on the leaf grid view.
     const auto& leafGridView = gridManager.grid().leafGridView();
 ```
 
 #### Set-up of the problem
-We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (scv) and
-subcontrolvolume faces (scvf) embedded in the elements of the grid partition.
+We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (`scv`) and
+subcontrolvolume faces (`scvf`) embedded in the elements of the grid partition.
 
 ```cpp
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
@@ -190,11 +190,11 @@ Convenience functions to perform a grid adaptation given a pre-calculated indica
         // In case of any adapted elements, the gridvariables and the pointsourcemap are updated.
         if (wasAdapted)
         {
-            // We overwrite the old solution with the new (resized & interpolated) one
+            // We overwrite the old solution with the new (resized & interpolated) one.
             xOld = x;
-            // We initialize the secondary variables to the new (and "new old") solution
+            // We initialize the secondary variables to the new (and "new old") solution.
             gridVariables->updateAfterGridAdaption(x);
-            // we update the point source map after adaption
+            // We update the point source map after adaption.
             problem->computePointSourceMap();
         }
     };
@@ -216,9 +216,8 @@ Afterwards, depending on the initial conditions, another grid adaptation using o
     adaptGridAndVariables(indicator);
 ```
 
-[[/codeblock]]
 #### Solving the problem
-We get some time loop parameters from the input file params.input
+We get some time loop parameters from the input file `params.input`
 
 ```cpp
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -227,7 +226,7 @@ We get some time loop parameters from the input file params.input
     const auto dt = getParam<Scalar>("TimeLoop.DtInitial");
 ```
 
-and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
+and initialize the VTK output. Each model has a predefined model-specific output with relevant parameters for that model.
 
 ```cpp
     using IOFields = GetPropType<TypeTag, Properties::IOFields>;
@@ -245,14 +244,14 @@ We instantiate the time loop
     timeLoop->setMaxTimeStepSize(maxDt);
 ```
 
-and set the assembler with the time loop because we have an instationary problem
+and set the assembler with the time loop because we have an instationary problem.
 
 ```cpp
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
 ```
 
-We set the linear solver and the non-linear solver
+We set the linear solver and the non-linear solver.
 
 ```cpp
     using LinearSolver = AMGBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
@@ -269,7 +268,7 @@ We start the time loop. In each time step before we start calculating a new solu
 ```cpp
     timeLoop->start(); do
     {
-        // We only want to refine/coarsen after first time step is finished, not before.
+        // We only want to refine/coarsen after the first time step is finished, not before.
         // The initial refinement was already done before the start of the time loop.
         // This means we only refine when the time is greater than 0. Note that we now also
         // have to update the assembler, since the sizes of the residual vector and jacobian matrix change.
@@ -290,20 +289,20 @@ We start the time loop. In each time step before we start calculating a new solu
         // We advance to the time loop to the next step.
         timeLoop->advanceTimeStep();
 
-        // We write vtk output for each time step
+        // We write vtk output for each time step.
         vtkWriter.write(timeLoop->time());
 
-        // We report statistics of this time step
+        // We report statistics of this time step.
         timeLoop->reportTimeStep();
 
-        // We set a new dt as suggested by the newton solver for the next time step
+        // We set a new dt as suggested by the newton solver for the next time step.
         timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
 
     } while (!timeLoop->finished());
 ```
 
 The following piece of code prints a final status report of the time loop
-before the program is terminated and we print he dumux end message
+before the program is terminated and we print he dumux end message.
 
 ```cpp
     timeLoop->finalize(leafGridView.comm());
diff --git a/examples/2pinfiltration/main.cc b/examples/2pinfiltration/main.cc
index 41db9c2cb599a2995180989e66f8e364bb166b11..bda04d93449004d31b906307a89f248ab4a80d9b 100644
--- a/examples/2pinfiltration/main.cc
+++ b/examples/2pinfiltration/main.cc
@@ -7,7 +7,7 @@
 // ## The file `main.cc`
 // [[content]]
 //
-// This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
+// This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method.
 // ### Included header files
 // [[details]] includes
 // [[codeblock]]
@@ -25,13 +25,13 @@
 #include <dumux/common/parameters.hh>
 #include <dumux/common/initialize.hh>
 
-//We include the linear solver to be used to solve the linear system and the nonlinear  Newton's method
+//We include the linear solver to be used to solve the linear system and the nonlinear Newton's method.
 #include <dumux/linear/istlsolvers.hh>
 #include <dumux/linear/linearsolvertraits.hh>
 #include <dumux/linear/linearalgebratraits.hh>
 #include <dumux/nonlinear/newtonsolver.hh>
 
-// Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a file that defines the different differentiation methods used to compute the derivatives of the residual
+// Further, we include the assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a header file that defines the different differentiation methods used to compute the derivatives of the residual.
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
 
@@ -40,14 +40,14 @@
 // The gridmanager constructs a grid from the information in the input or grid file.
 #include <dumux/io/grid/gridmanager_alu.hh>
 
-//We include several files which are needed for the adaptive grid
+//We include several files which are needed for the adaptive grid.
 #include <dumux/adaptive/adapt.hh>
 #include <dumux/adaptive/markelements.hh>
 #include <dumux/adaptive/initializationindicator.hh>
 #include <dumux/porousmediumflow/2p/griddatatransfer.hh>
 #include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
 
-// Finally, we include the properties which configure the simulation
+// Finally, we include the properties which configure the simulation.
 #include "properties.hh"
 // [[/details]]
 //
@@ -81,13 +81,13 @@ int main(int argc, char** argv) try
 
     // The instationary non-linear problem is run on this grid.
     //
-    // we compute on the leaf grid view
+    // We compute on the leaf grid view.
     const auto& leafGridView = gridManager.grid().leafGridView();
     // [[/codeblock]]
 
     // #### Set-up of the problem
-    // We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (scv) and
-    // subcontrolvolume faces (scvf) embedded in the elements of the grid partition.
+    // We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (`scv`) and
+    // subcontrolvolume faces (`scvf`) embedded in the elements of the grid partition.
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
     auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
 
@@ -140,11 +140,11 @@ int main(int argc, char** argv) try
         // In case of any adapted elements, the gridvariables and the pointsourcemap are updated.
         if (wasAdapted)
         {
-            // We overwrite the old solution with the new (resized & interpolated) one
+            // We overwrite the old solution with the new (resized & interpolated) one.
             xOld = x;
-            // We initialize the secondary variables to the new (and "new old") solution
+            // We initialize the secondary variables to the new (and "new old") solution.
             gridVariables->updateAfterGridAdaption(x);
-            // we update the point source map after adaption
+            // We update the point source map after adaption.
             problem->computePointSourceMap();
         }
     };
@@ -153,6 +153,7 @@ int main(int argc, char** argv) try
     // We do initial refinements around sources/BCs, for which we use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh`.
     // Afterwards, depending on the initial conditions, another grid adaptation using our `indicator` above might be necessary, which uses
     // `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
+    // [[codeblock]]
     GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
     const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
     for (std::size_t i = 0; i < maxLevel; ++i)
@@ -166,13 +167,13 @@ int main(int argc, char** argv) try
 
     // #### Solving the problem
 
-    // We get some time loop parameters from the input file params.input
+    // We get some time loop parameters from the input file `params.input`
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
     const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
     const auto dt = getParam<Scalar>("TimeLoop.DtInitial");
 
-    // and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
+    // and initialize the VTK output. Each model has a predefined model-specific output with relevant parameters for that model.
     using IOFields = GetPropType<TypeTag, Properties::IOFields>;
     VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
     using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
@@ -184,11 +185,11 @@ int main(int argc, char** argv) try
     auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 
-    // and set the assembler with the time loop because we have an instationary problem
+    // and set the assembler with the time loop because we have an instationary problem.
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
 
-    // We set the linear solver and the non-linear solver
+    // We set the linear solver and the non-linear solver.
     using LinearSolver = AMGBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
                                                LinearAlgebraTraitsFromAssembler<Assembler>>;
     auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
@@ -201,7 +202,7 @@ int main(int argc, char** argv) try
     // [[codeblock]]
     timeLoop->start(); do
     {
-        // We only want to refine/coarsen after first time step is finished, not before.
+        // We only want to refine/coarsen after the first time step is finished, not before.
         // The initial refinement was already done before the start of the time loop.
         // This means we only refine when the time is greater than 0. Note that we now also
         // have to update the assembler, since the sizes of the residual vector and jacobian matrix change.
@@ -222,20 +223,20 @@ int main(int argc, char** argv) try
         // We advance to the time loop to the next step.
         timeLoop->advanceTimeStep();
 
-        // We write vtk output for each time step
+        // We write vtk output for each time step.
         vtkWriter.write(timeLoop->time());
 
-        // We report statistics of this time step
+        // We report statistics of this time step.
         timeLoop->reportTimeStep();
 
-        // We set a new dt as suggested by the newton solver for the next time step
+        // We set a new dt as suggested by the newton solver for the next time step.
         timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
 
     } while (!timeLoop->finished());
     // [[/codeblock]]
 
     // The following piece of code prints a final status report of the time loop
-    //  before the program is terminated and we print he dumux end message
+    //  before the program is terminated and we print he dumux end message.
     // [[codeblock]]
     timeLoop->finalize(leafGridView.comm());
 
diff --git a/examples/2pinfiltration/problem.hh b/examples/2pinfiltration/problem.hh
index 3c47952304ec385382119b246729f323999026e1..d90a16ea8ced75736b21b5630cf18f73c9b1ec11 100644
--- a/examples/2pinfiltration/problem.hh
+++ b/examples/2pinfiltration/problem.hh
@@ -117,8 +117,8 @@ public:
 
     // #### Neumann boundaries
     // In our case, we need to specify mass fluxes for our two liquid phases.
-    // Negative sign means influx and the unit of the boundary flux is $`kg/(m^2 s)`$.
-    // On the inlet area, we set a DNAPL influx of $`0.04 kg/(m^2 s)`$. On all other
+    // Negative sign means influx and the unit of the boundary flux is $\mathrm{kg}/\mathrm{(m^2 s)}$.
+    // On the inlet area, we set a DNAPL influx of $0.04 \, \mathrm{kg}/\mathrm{(m^2 s)}$. On all other
     // Neumann boundaries, the boundary flux is zero.
     // [[codeblock]]
     NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
@@ -154,7 +154,7 @@ public:
     // Point sources are added by pushing them into the vector `pointSources`.
     // The `PointSource` constructor takes two arguments.
     // The first argument is a coordinate array containing the position in space,
-    // the second argument is an array of source value for each equation (in units of $`kg/s`$).
+    // the second argument is an array of source values for each equation (in units of $`\mathrm{kg}/\mathrm{s}`$).
     // Recall that the first equation is the water phase mass balance
     // and the second equation is the DNAPL phase mass balance.
     void addPointSources(std::vector<PointSource>& pointSources) const
diff --git a/examples/2pinfiltration/properties.hh b/examples/2pinfiltration/properties.hh
index 00d2f93d646c50997fcc5d9f2178d434bbb8609e..2301b779296a1d64419d6490ec842016964f2e1a 100644
--- a/examples/2pinfiltration/properties.hh
+++ b/examples/2pinfiltration/properties.hh
@@ -72,7 +72,7 @@ struct SpatialParams<TypeTag, TTag::PointSourceExample>
 // The `Grid` property tells the
 // simulator to use ALUGrid - an unstructured grid manager - here
 // configured for grid and coordinate dimensions `2`,
-// hexahedral element types (`Dune::cube`) and non-conforming refinement mode.
+// hexahedral element types (`Dune::cube`) and non-conforming refinement mode (`Dune::nonconforming`).
 // `Dune::ALUGrid` is declared in the included header `dune/alugrid/grid.hh`
 // from the Dune module `dune-alugrid`.
 template<class TypeTag>
diff --git a/examples/biomineralization/doc/mainfile.md b/examples/biomineralization/doc/mainfile.md
index a4bf1bd210a625a92d78007a3c1b0ce14942ca28..9998eacf9eee776cba524f41142c261b9d90a224 100644
--- a/examples/biomineralization/doc/mainfile.md
+++ b/examples/biomineralization/doc/mainfile.md
@@ -6,7 +6,7 @@
 
 # Main file
 
-The main file features a time loo with check points to assure the correct timing of the injections matching an experimental setup (see [@Hommel2016], Section 4.2 under the name _column experiment D1_)
+The main file features a time loop with check points to assure the correct timing of the injections matching an experimental setup (see [@Hommel2016], Section 4.2 under the name _column experiment D1_).
 The timing of the injections is stored in an additional input file `injections_checkpoints.dat`,
 which is read by `main.cc` to set the check points in the time loop.
 The number of check points reached is counted and passed on to `problem.hh`, so that `problem.hh` can determine the appropriate injection type.
@@ -167,9 +167,9 @@ We initialize the vtk output module. Each model has a predefined model specific
     using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
     vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
     GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter); //!< Add model specific output fields
-    //we add permeability as a specific output
+    // We add permeability as a specific output
     vtkWriter.addField(problem->getPermeability(), "Permeability");
-    //We update the output fields before write
+    // We update the output fields before write
     problem->updateVtkOutput(x);
     vtkWriter.write(0.0);
 ```
@@ -178,24 +178,24 @@ We instantiate the time loop and define an episode index to keep track of the th
 
 ```cpp
     int episodeIdx = 0;
-    //We  set the initial episodeIdx in the problem to zero
+    // We set the initial episodeIdx in the problem to zero
     problem->setEpisodeIdx(episodeIdx);
-    //We  set the time loop with episodes (check points)
+    // We set the time loop with episodes (check points)
     auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 ```
 
 ### Step 3 Setting episodes specified from file
 In this example we want to specify the injected reactant solutions at certain times.
-This is specified in an external file injections_checkpoints.dat.
+This is specified in an external file `injections_checkpoints.dat`.
 Within the following code block, the parameter file read and the time for the injection are extracted.
 Based on that, the checkpoints for the simulation are set.
 
 ```cpp
-    //We use a time loop with episodes if an injection parameter file is specified in the input
+    // We use a time loop with episodes if an injection parameter file is specified in the input
     if (hasParam("Injection.NumInjections"))
     {
-        //We first read the times of the checkpoints from the injection file
+        // We first read the times of the checkpoints from the injection file
         const auto injectionCheckPoints = readFileToContainer<std::vector<double>>("injection_checkpoints.dat");
 
         // We set the time loop with episodes of various lengths as specified in the injection file
@@ -214,7 +214,7 @@ Based on that, the checkpoints for the simulation are set.
     }
 ```
 
-### Step 4 solving the instationary problem
+### Step 4 Solving the instationary problem
 We create and initialize the assembler with a time loop for the transient problem.
 Within the time loop, we will use this assembler in each time step to assemble the linear system.
 Additionally the linear and non-linear solvers are set
@@ -223,13 +223,13 @@ Additionally the linear and non-linear solvers are set
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
 
-    //We set the linear solver
+    // We set the linear solver
     using LinearSolver = ILUBiCGSTABIstlSolver<
         LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>
     >;
     auto linearSolver = std::make_shared<LinearSolver>();
 
-    //We set the non-linear solver
+    // We set the non-linear solver
     using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 ```
@@ -270,7 +270,7 @@ and the time step sizes used is printed to the terminal.
         // We report statistics of this time step
         timeLoop->reportTimeStep();
 
-        //If episodes/check points are used, we count episodes and update the episode indices in the main file and the problem.
+        // If episodes/check points are used, we count episodes and update the episode indices in the main file and the problem.
         if (hasParam("Injection.NumInjections") || hasParam("TimeLoop.EpisodeLength"))
         {
             if (timeLoop->isCheckPoint())
diff --git a/examples/biomineralization/doc/mainfile_intro.md b/examples/biomineralization/doc/mainfile_intro.md
index c1435eb8998641ec2ee730ae02646974ef5ed92a..1a67f6f312c0fc443ac49d5a5505a3266b3f854d 100644
--- a/examples/biomineralization/doc/mainfile_intro.md
+++ b/examples/biomineralization/doc/mainfile_intro.md
@@ -1,6 +1,6 @@
 # Main file
 
-The main file features a time loo with check points to assure the correct timing of the injections matching an experimental setup (see [@Hommel2016], Section 4.2 under the name _column experiment D1_)
+The main file features a time loop with check points to assure the correct timing of the injections matching an experimental setup (see [@Hommel2016], Section 4.2 under the name _column experiment D1_).
 The timing of the injections is stored in an additional input file `injections_checkpoints.dat`,
 which is read by `main.cc` to set the check points in the time loop.
 The number of check points reached is counted and passed on to `problem.hh`, so that `problem.hh` can determine the appropriate injection type.
diff --git a/examples/biomineralization/doc/modelconcept.md b/examples/biomineralization/doc/modelconcept.md
index d31c4b6f218b54e279c1ea5d5e8420a976922dc4..48b67319f2302b6f94afba916fafb6aed1597a60 100644
--- a/examples/biomineralization/doc/modelconcept.md
+++ b/examples/biomineralization/doc/modelconcept.md
@@ -13,7 +13,7 @@ A mass transfer may occur between both fluid phases by the mutual dissolution of
 the aqueous phase and the two *solid* phases, biofilm and calcite, denoted by subscripts (f) and (c) respectively, by attachment or
 detachment of biomass and precipitation or dissolution of calcite.
 
-The mobile components, denoted by superscripts k, are water (w), dissolved inorganic carbon (C<sub>tot</sub>),
+The mobile components, denoted by superscripts $`\kappa`$, are water (w), dissolved inorganic carbon (C<sub>tot</sub>),
 sodium (Na), chloride (Cl), calcium (Ca), urea (u), substrate (s), oxygen (O<sub>2</sub>), and suspended biomass (b).
 A substrate is the carbon and energy source of the bacteria and O<sub>2</sub> is the electron acceptor.
 
@@ -52,8 +52,8 @@ only a storage and source term since they are immobile:
 ```
 
 Here, $`\phi_\lambda`$ and $`\rho_\lambda`$ are volume fraction and mass density of
-the solid phase $`\lambda`$, and $`q^\lambda`$ is the source term of phase $`\lambda`$ due to biochemical reactions.
-The sources and sinks due to reactions $`q^\kappa$ and $q^\lambda`$ are specific to the components and are discussed in details in the subsequent section.
+the solid phase $`\lambda`$, and $`q^{\lambda}`$ is the source term of phase $`\lambda`$ due to biochemical reactions.
+The sources and sinks due to reactions $`q^{\kappa}$ and $q^{\lambda}`$ are specific to the components and are discussed in details in the subsequent section.
 
 ## Component-specific reactive source and sink terms
 
diff --git a/examples/biomineralization/doc/modelconcept_intro.md b/examples/biomineralization/doc/modelconcept_intro.md
index fb953e4992b6361b8bc7150a7aafb9e7d37aa7c2..af9a71edf6807422f30de13fc5f651a105d1c11d 100644
--- a/examples/biomineralization/doc/modelconcept_intro.md
+++ b/examples/biomineralization/doc/modelconcept_intro.md
@@ -7,7 +7,7 @@ A mass transfer may occur between both fluid phases by the mutual dissolution of
 the aqueous phase and the two *solid* phases, biofilm and calcite, denoted by subscripts (f) and (c) respectively, by attachment or
 detachment of biomass and precipitation or dissolution of calcite.
 
-The mobile components, denoted by superscripts k, are water (w), dissolved inorganic carbon (C<sub>tot</sub>),
+The mobile components, denoted by superscripts $`\kappa`$, are water (w), dissolved inorganic carbon (C<sub>tot</sub>),
 sodium (Na), chloride (Cl), calcium (Ca), urea (u), substrate (s), oxygen (O<sub>2</sub>), and suspended biomass (b).
 A substrate is the carbon and energy source of the bacteria and O<sub>2</sub> is the electron acceptor.
 
@@ -46,8 +46,8 @@ only a storage and source term since they are immobile:
 ```
 
 Here, $`\phi_\lambda`$ and $`\rho_\lambda`$ are volume fraction and mass density of
-the solid phase $`\lambda`$, and $`q^\lambda`$ is the source term of phase $`\lambda`$ due to biochemical reactions.
-The sources and sinks due to reactions $`q^\kappa$ and $q^\lambda`$ are specific to the components and are discussed in details in the subsequent section.
+the solid phase $`\lambda`$, and $`q^{\lambda}`$ is the source term of phase $`\lambda`$ due to biochemical reactions.
+The sources and sinks due to reactions $`q^{\kappa}$ and $q^{\lambda}`$ are specific to the components and are discussed in details in the subsequent section.
 
 ## Component-specific reactive source and sink terms
 
diff --git a/examples/biomineralization/doc/setup.md b/examples/biomineralization/doc/setup.md
index 52f5a8ec59b4025d8882d0bcabc1e5124f05b53b..2e53bc5bbed6c05e99c8fa04332f87a95ec68bf1 100644
--- a/examples/biomineralization/doc/setup.md
+++ b/examples/biomineralization/doc/setup.md
@@ -808,7 +808,7 @@ for which we then specialize `properties` to the needs of this setup.
 <details>
 
 This type tag specializes most of the `properties` required for two phase flow with
-multiple components including mineralisation simulations (2pncmin) in DuMuX
+multiple components including mineralisation simulations (2pncmin) in DuMu<sup>x</sup>
 We will use this in the following to inherit the respective properties and
 subsequently specialize those `properties` for our `TypeTag`, which we want to
 modify or for which no meaningful default can be set.
diff --git a/examples/biomineralization/main.cc b/examples/biomineralization/main.cc
index fc32029128cbb3d0d9420b61d16706badc4f0800..efbb88c9ccc6e7921f44f152c513c02dbf65deb7 100644
--- a/examples/biomineralization/main.cc
+++ b/examples/biomineralization/main.cc
@@ -116,34 +116,34 @@ int main(int argc, char** argv) try
     using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
     vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
     GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter); //!< Add model specific output fields
-    //we add permeability as a specific output
+    // We add permeability as a specific output
     vtkWriter.addField(problem->getPermeability(), "Permeability");
-    //We update the output fields before write
+    // We update the output fields before write
     problem->updateVtkOutput(x);
     vtkWriter.write(0.0);
     // [[/codeblock]]
 
-    //We instantiate the time loop and define an episode index to keep track of the the actual episode number
+    // We instantiate the time loop and define an episode index to keep track of the the actual episode number
     // [[codeblock]]
     int episodeIdx = 0;
-    //We  set the initial episodeIdx in the problem to zero
+    // We set the initial episodeIdx in the problem to zero
     problem->setEpisodeIdx(episodeIdx);
-    //We  set the time loop with episodes (check points)
+    // We set the time loop with episodes (check points)
     auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
     // [[/codeblock]]
 
     // ### Step 3 Setting episodes specified from file
     // In this example we want to specify the injected reactant solutions at certain times.
-    // This is specified in an external file injections_checkpoints.dat.
+    // This is specified in an external file `injections_checkpoints.dat`.
     // Within the following code block, the parameter file read and the time for the injection are extracted.
     // Based on that, the checkpoints for the simulation are set.
 
     // [[codeblock]]
-    //We use a time loop with episodes if an injection parameter file is specified in the input
+    // We use a time loop with episodes if an injection parameter file is specified in the input
     if (hasParam("Injection.NumInjections"))
     {
-        //We first read the times of the checkpoints from the injection file
+        // We first read the times of the checkpoints from the injection file
         const auto injectionCheckPoints = readFileToContainer<std::vector<double>>("injection_checkpoints.dat");
 
         // We set the time loop with episodes of various lengths as specified in the injection file
@@ -162,7 +162,7 @@ int main(int argc, char** argv) try
     }
     // [[/codeblock]]
 
-    // ### Step 4 solving the instationary problem
+    // ### Step 4 Solving the instationary problem
 
     // We create and initialize the assembler with a time loop for the transient problem.
     // Within the time loop, we will use this assembler in each time step to assemble the linear system.
@@ -171,13 +171,13 @@ int main(int argc, char** argv) try
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
 
-    //We set the linear solver
+    // We set the linear solver
     using LinearSolver = ILUBiCGSTABIstlSolver<
         LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>
     >;
     auto linearSolver = std::make_shared<LinearSolver>();
 
-    //We set the non-linear solver
+    // We set the non-linear solver
     using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
     // [[/codeblock]]
@@ -218,7 +218,7 @@ int main(int argc, char** argv) try
         // We report statistics of this time step
         timeLoop->reportTimeStep();
 
-        //If episodes/check points are used, we count episodes and update the episode indices in the main file and the problem.
+        // If episodes/check points are used, we count episodes and update the episode indices in the main file and the problem.
         if (hasParam("Injection.NumInjections") || hasParam("TimeLoop.EpisodeLength"))
         {
             if (timeLoop->isCheckPoint())
diff --git a/examples/biomineralization/properties.hh b/examples/biomineralization/properties.hh
index 323eead9e009f0b758e9380d336661afc86fb178..604435af44471045276c148c5cfea2600ffbcae9 100644
--- a/examples/biomineralization/properties.hh
+++ b/examples/biomineralization/properties.hh
@@ -19,7 +19,7 @@
 // <details>
 //
 // This type tag specializes most of the `properties` required for two phase flow with
-// multiple components including mineralisation simulations (2pncmin) in DuMuX
+// multiple components including mineralisation simulations (2pncmin) in DuMu<sup>x</sup>
 // We will use this in the following to inherit the respective properties and
 // subsequently specialize those `properties` for our `TypeTag`, which we want to
 // modify or for which no meaningful default can be set.
diff --git a/examples/embedded_network_1d3d/README.md b/examples/embedded_network_1d3d/README.md
index 599eb3565c134e4f65fcd8c161abc2d4761f927c..e4ac11f676c848bc534a78ec24214e1920385dbe 100644
--- a/examples/embedded_network_1d3d/README.md
+++ b/examples/embedded_network_1d3d/README.md
@@ -2,7 +2,7 @@
 
 # Embedded network 1D-3D model (tissue perfusion)
 
-We solve a tracer transport problem in a domain that consist of a porous medium block
+We solve a tracer transport problem in a domain that consists of a porous medium block
 that has an embedded transport network. In this example the porous medium is brain tissue
 and the embedded network is the vasculature (blood vessels).
 
@@ -12,8 +12,8 @@ __Table of contents__. This description is structured as follows:
 
 ## Problem set-up
 In this example we simulate clearance of a substance present in the tissue through the blood. The tissue
-cube is assigned no-flux/symmetry boundary conditions assuming that identical cubes are mirrored on all sides.
-Therefore, the tracer has to cross the vessel wall into the network (vessel lumen). It then gets transported
+cube is assigned no-flux/symmetry boundary conditions, assuming that identical cubes are mirrored on all sides.
+Therefore, the tracer has to cross the vessel wall into the network (vessel lumen). It is then transported
 in the blood stream by advection and diffusion. In the network, the tracer mole fraction is zero at the inlet
 and at the outlet the mole fraction gradient is zero. Thus, the tracer is transported out of the domain by advection only.
 VTK output is written in every time step, and the total tracer concentration in the tissue is written into a text file
@@ -23,7 +23,7 @@ along the simulation.
 The domain consists of a small blood vessel network embedded
 in a porous tissue (pore space is the pore space of the extra-cellular matrix).
 The network is extracted based on the mouse cortical brain data from [Blinder (2013)](https://doi.org/10.1038/nn.3426).
-With the boundary condition estimated by [Schmid (2017a)](https://doi.org/10.1371/journal.pcbi.1005392)
+With the boundary conditions estimated by [Schmid (2017a)](https://doi.org/10.1371/journal.pcbi.1005392)
 and the pressure and network data available from [Schmid (2017b)](https://doi.org/10.5281/zenodo.269650), blood flow
 simulations in DuMu<sup>x</sup> give a pressure field in the entire network. A small part (200µm)³ has been extracted
 for this example. The example also contains a blood flow solver which uses the pressure boundary data to compute
@@ -56,7 +56,7 @@ $`Q_\mathrm{B} := A_\mathrm{B}v_\mathrm{B}`$ is a given blood flow field transpo
 Furthermore, isothermal conditions with a homogeneous temperature distribution of constant $`T=37^\circ C`$ are assumed.
 The 1D network PDE is formulated in terms of the local axial coordinate $`s`$.
 
-For the coupling, we use a formulation with line source and a perimeter integral operator, that is $`\Phi_\Lambda`$ is chosen to be a delta distribution centered on vessel centerline, and $`x_\mathrm{B}`$, which is assumed to be constant on each cross-section (well-mixed),
+For the coupling, we use a formulation with line source and a perimeter integral operator, denoted  $`\Phi_\Lambda`$, which is chosen to be a delta distribution centered on vessel centerline, and $`x_\mathrm{B}`$, which is assumed to be constant on each cross-section (well-mixed),
 and every point on the surface is assumed to be uniquely identified with a position $`s`$ in the network such that we can extend
 $`x_\mathrm{B}`$ with $`\Pi`$ to the surface. The integral is evaluated numerically. Each integration point is represented in
 the code as a point source. The point source values (the integrand) are implemented in the problem class function `pointSource(...)`.
diff --git a/examples/embedded_network_1d3d/doc/_intro.md b/examples/embedded_network_1d3d/doc/_intro.md
index 32aaf2b9037763edb8f522524d208843793414d5..a93e21a735d395f1fdbd425e3c68ca243a61a87a 100644
--- a/examples/embedded_network_1d3d/doc/_intro.md
+++ b/examples/embedded_network_1d3d/doc/_intro.md
@@ -1,6 +1,6 @@
 # Embedded network 1D-3D model (tissue perfusion)
 
-We solve a tracer transport problem in a domain that consist of a porous medium block
+We solve a tracer transport problem in a domain that consists of a porous medium block
 that has an embedded transport network. In this example the porous medium is brain tissue
 and the embedded network is the vasculature (blood vessels).
 
@@ -10,8 +10,8 @@ __Table of contents__. This description is structured as follows:
 
 ## Problem set-up
 In this example we simulate clearance of a substance present in the tissue through the blood. The tissue
-cube is assigned no-flux/symmetry boundary conditions assuming that identical cubes are mirrored on all sides.
-Therefore, the tracer has to cross the vessel wall into the network (vessel lumen). It then gets transported
+cube is assigned no-flux/symmetry boundary conditions, assuming that identical cubes are mirrored on all sides.
+Therefore, the tracer has to cross the vessel wall into the network (vessel lumen). It is then transported
 in the blood stream by advection and diffusion. In the network, the tracer mole fraction is zero at the inlet
 and at the outlet the mole fraction gradient is zero. Thus, the tracer is transported out of the domain by advection only.
 VTK output is written in every time step, and the total tracer concentration in the tissue is written into a text file
@@ -21,7 +21,7 @@ along the simulation.
 The domain consists of a small blood vessel network embedded
 in a porous tissue (pore space is the pore space of the extra-cellular matrix).
 The network is extracted based on the mouse cortical brain data from [Blinder (2013)](https://doi.org/10.1038/nn.3426).
-With the boundary condition estimated by [Schmid (2017a)](https://doi.org/10.1371/journal.pcbi.1005392)
+With the boundary conditions estimated by [Schmid (2017a)](https://doi.org/10.1371/journal.pcbi.1005392)
 and the pressure and network data available from [Schmid (2017b)](https://doi.org/10.5281/zenodo.269650), blood flow
 simulations in DuMu<sup>x</sup> give a pressure field in the entire network. A small part (200µm)³ has been extracted
 for this example. The example also contains a blood flow solver which uses the pressure boundary data to compute
@@ -54,7 +54,7 @@ $`Q_\mathrm{B} := A_\mathrm{B}v_\mathrm{B}`$ is a given blood flow field transpo
 Furthermore, isothermal conditions with a homogeneous temperature distribution of constant $`T=37^\circ C`$ are assumed.
 The 1D network PDE is formulated in terms of the local axial coordinate $`s`$.
 
-For the coupling, we use a formulation with line source and a perimeter integral operator, that is $`\Phi_\Lambda`$ is chosen to be a delta distribution centered on vessel centerline, and $`x_\mathrm{B}`$, which is assumed to be constant on each cross-section (well-mixed),
+For the coupling, we use a formulation with line source and a perimeter integral operator, denoted  $`\Phi_\Lambda`$, which is chosen to be a delta distribution centered on vessel centerline, and $`x_\mathrm{B}`$, which is assumed to be constant on each cross-section (well-mixed),
 and every point on the surface is assumed to be uniquely identified with a position $`s`$ in the network such that we can extend
 $`x_\mathrm{B}`$ with $`\Pi`$ to the surface. The integral is evaluated numerically. Each integration point is represented in
 the code as a point source. The point source values (the integrand) are implemented in the problem class function `pointSource(...)`.
diff --git a/examples/embedded_network_1d3d/doc/main.md b/examples/embedded_network_1d3d/doc/main.md
index bb4b5683a7998668e032a39b2f1e2074c68bc428..f9e67a16662799d3f1bc3ba289fcaf700b412f6d 100644
--- a/examples/embedded_network_1d3d/doc/main.md
+++ b/examples/embedded_network_1d3d/doc/main.md
@@ -11,7 +11,7 @@
 This file contains the `main` function and implements the main program
 flow. We initialize the simulation framework, initialize parameters,
 create the grids (using parameters from the configuration file `params.input`),
-create vector to store the primary and secondary variables, construct an
+create vectors to store the primary and secondary variables, construct an
 assembler that can assembly the the residual (discrete equations) and
 the system matrix (Jacobian of the residual), and create a linear solver
 that solves the linear system arising in each time step. The time loop
diff --git a/examples/embedded_network_1d3d/doc/spatialparams.md b/examples/embedded_network_1d3d/doc/spatialparams.md
index f4379d7773cb4d6ca93234cd55086ea718176111..3b1cf1d1353681cd2aca90a547c53521ae75df2c 100644
--- a/examples/embedded_network_1d3d/doc/spatialparams.md
+++ b/examples/embedded_network_1d3d/doc/spatialparams.md
@@ -178,7 +178,7 @@ public:
 
 
 Some interface functions that can be used, for example, in the problem class, to
-access spatial parameter. The interface `radius` is required by the 1D-3D coupling
+access spatial parameters. The interface `radius` is required by the 1D-3D coupling
 manager and is used to create the integration points for the coupling term
 
 
diff --git a/examples/embedded_network_1d3d/main.cc b/examples/embedded_network_1d3d/main.cc
index 97c0ab3653bc73a394d795b5878af5757659f1d7..d9519f69308254d89c7a47d06601cf7d45fba783 100644
--- a/examples/embedded_network_1d3d/main.cc
+++ b/examples/embedded_network_1d3d/main.cc
@@ -9,7 +9,7 @@
 // This file contains the `main` function and implements the main program
 // flow. We initialize the simulation framework, initialize parameters,
 // create the grids (using parameters from the configuration file `params.input`),
-// create vector to store the primary and secondary variables, construct an
+// create vectors to store the primary and secondary variables, construct an
 // assembler that can assembly the the residual (discrete equations) and
 // the system matrix (Jacobian of the residual), and create a linear solver
 // that solves the linear system arising in each time step. The time loop
diff --git a/examples/embedded_network_1d3d/spatialparams.hh b/examples/embedded_network_1d3d/spatialparams.hh
index 4d36e0ae01e1125e078a9eac8ddcc6fe95429fee..dfb7a8c2bffb48555445437aa3c3e769c560a57f 100644
--- a/examples/embedded_network_1d3d/spatialparams.hh
+++ b/examples/embedded_network_1d3d/spatialparams.hh
@@ -157,7 +157,7 @@ public:
     // [[/codeblock]]
     //
     // Some interface functions that can be used, for example, in the problem class, to
-    // access spatial parameter. The interface `radius` is required by the 1D-3D coupling
+    // access spatial parameters. The interface `radius` is required by the 1D-3D coupling
     // manager and is used to create the integration points for the coupling term
     //
     // [[codeblock]]
diff --git a/examples/freeflowchannel/README.md b/examples/freeflowchannel/README.md
index db1040f77ffb1a3dd5af6ba9b336073cfcd04762..29b7b0013a02ce4f3ddacb7c0ff795b8d66daa37 100644
--- a/examples/freeflowchannel/README.md
+++ b/examples/freeflowchannel/README.md
@@ -41,7 +41,7 @@ have a look at the DuMu<sup>x</sup> [documentation](https://dumux.org/docs/doxyg
 ## Problem set-up
 This example considers stationary flow of a fluid between two parallel solid plates in two dimensions.
 Flow is enforced from left to right by prescribing an inflow velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$
-on the left boundary, while a fixed pressure of $`p = 1.1 \text{bar}`$ and a zero velocity gradient
+on the left boundary, while a fixed pressure of $`p = 1.1 \, \text{bar}`$ and a zero velocity gradient
 in x-direction are prescribed on the right boundary. On the top and bottom boundaries, no-slip
 conditions are applied, which cause a parabolic velocity profile to develop along the channel.
 Take a look at Figure 1 for an illustration of the domain and the boundary conditions.
@@ -69,7 +69,7 @@ Take a look at Figure 1 for an illustration of the domain and the boundary condi
 
 ## Compile-time settings (`properties.hh`)
 
-In this file, the type tag used for this simulation is defined,
+In this file, the type tag used for this simulation (`TTag::ChannelExample`) is defined,
 for which we then specialize properties (compile time options) to the needs of the desired setup.
 
 
@@ -81,7 +81,7 @@ for which we then specialize properties (compile time options) to the needs of t
 <details><summary> Click to show includes</summary>
 
 The `NavierStokes` type tag specializes most of the properties required for Navier-
-Stokes single-phase flow simulations in DuMuX. We will use this in the following to inherit the
+Stokes single-phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
 respective properties and subsequently specialize those properties for our
 type tag, which we want to modify or for which no meaningful default can be set.
 
@@ -132,7 +132,7 @@ Navier-Stokes equations. This can be achieved at runtime by setting the paramete
 to see how this is done in this example.
 
 ```cpp
-// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
+// We enter the namespace `Dumux::Properties` in order to import the entire `Dumux` namespace for general use:
 namespace Dumux::Properties {
 
 // declaration of the `ChannelExample` type tag for the single-phase flow problem
@@ -151,11 +151,11 @@ default can be set, are specialized for our type tag `ChannelExample`.
 template<class TypeTag>
 struct Grid<TypeTag, TTag::ChannelExample> { using type = Dune::YaspGrid<2>; };
 
-// This sets our problem class (see problem.hh) containing initial and boundary conditions.
+// This sets our problem type (see `problem.hh`) containing the initial and boundary conditions.
 template<class TypeTag>
 struct Problem<TypeTag, TTag::ChannelExample> { using type = Dumux::ChannelExampleProblem<TypeTag> ; };
 
-// This sets the fluid system to be used. Here, we use a liquid with constant properties as fluid phase.
+// This sets the fluid system type to be used. Here, we use a liquid with constant properties as the fluid phase.
 template<class TypeTag>
 struct FluidSystem<TypeTag, TTag::ChannelExample>
 {
@@ -168,7 +168,7 @@ We also set some properties related to memory management
 throughout the simulation.
 <details><summary> Click to show caching properties</summary>
 
-In Dumux, one has the option to activate/deactivate the grid-wide caching of
+In DuMu<sup>x</sup>, one has the option to activate/deactivate the grid-wide caching of
 geometries and variables. If active, the CPU time can be significantly reduced
 as less dynamic memory allocation procedures are necessary. Per default, grid-wide
 caching is disabled to ensure minimal memory requirements, however, in this example we
@@ -221,7 +221,7 @@ Include the `NavierStokesBoundaryTypes` class which specifies the boundary types
 ```
 
 ### The problem class
-We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+We enter the problem class `ChannelExampleProblem` where all necessary boundary conditions and initial conditions are set for our simulation.
 As we are solving a problem related to free flow, we inherit from the base class `NavierStokesStaggeredProblem`.
 
 ```cpp
@@ -363,7 +363,7 @@ Finally, private variables are declared:
 
 ### Included header files
 <details><summary> Click to show includes</summary>
-These are DUNE helper classes related to parallel computations and file I/O
+These are DUNE helper classes related to parallel computations and file I/O.
 
 ```cpp
 #include <dune/common/parallel/mpihelper.hh>
@@ -442,7 +442,7 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
 ```
 
-We define a convenience alias for the type tag of the problem. The type
+We define a convenience alias for the type tag of the problem (`Properties::TTag::ChannelExample`). The type
 tag contains all the properties that are needed to define the model and the problem
 setup. Throughout the main file, we will obtain types defined for this type tag
 using the property system, i.e. with `GetPropType`.
@@ -466,7 +466,7 @@ of the corners of the grid and the number of cells to be used to discretize each
 
 #### Step 2: Setting up and solving the problem
 First, a finite volume grid geometry is constructed from the grid that was created above.
-This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
+This builds the sub-control volumes (`scv`) and sub-control volume faces (`scvf`) for each element
 of the grid partition.
 
 ```cpp
@@ -474,14 +474,14 @@ of the grid partition.
     auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
 ```
 
-We now instantiate the problem, in which we define the boundary and initial conditions.
+We now instantiate the `problem`, in which we define the boundary and initial conditions.
 
 ```cpp
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     auto problem = std::make_shared<Problem>(gridGeometry);
 ```
 
-We set a solution vector which consist of two parts: one part (indexed by `cellCenterIdx`)
+We set a solution vector `x` which consist of two parts: one part (indexed by `cellCenterIdx`)
 is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
 (indexed by `faceIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
 We initialize the solution vector by what was defined as the initial solution of the the problem.
@@ -497,12 +497,12 @@ We initialize the solution vector by what was defined as the initial solution of
 The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
 
 ```cpp
-    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using GridVariables =GetPropType<TypeTag, Properties::GridVariables>;
     auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
     gridVariables->init(x);
 ```
 
-We then initialize the predefined model-specific output vtk output.
+We then initialize the predefined model-specific output VTK output.
 
 ```cpp
     using IOFields = GetPropType<TypeTag, Properties::IOFields>;
@@ -513,13 +513,13 @@ We then initialize the predefined model-specific output vtk output.
 
 <details><summary> Click to show calculation of surface fluxes</summary>
 We set up two surfaces over which fluxes are calculated.
-We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain.
-The first surface (added by the first call of addSurface) shall be placed at the middle of the channel.
+We determine the extend $`[xMin,xMax] \times [yMin,yMax]`$ of the physical domain.
+The first surface (added by the first call of `addSurface`) shall be placed at the middle of the channel.
 If we have an odd number of cells in x-direction, there would not be any cell faces
 at the position of the surface (which is required for the flux calculation).
 In this case, we add half a cell-width to the x-position in order to make sure that
 the cell faces lie on the surface. This assumes a regular cartesian grid.
-The second surface (second call of addSurface) is placed at the outlet of the channel.
+The second surface (second call of `addSurface`) is placed at the outlet of the channel.
 
 ```cpp
     FluxOverSurface<GridVariables,
@@ -556,7 +556,7 @@ The second surface (second call of addSurface) is placed at the outlet of the ch
 ```
 
 </details>
-We create and initialize the assembler for the stationary problem.
+We create and initialize the `assembler` for the stationary problem.
 This is where the Jacobian matrix for the Newton solver is assembled.
 
 ```cpp
diff --git a/examples/freeflowchannel/doc/_intro.md b/examples/freeflowchannel/doc/_intro.md
index d3325394f4043140dc6d66ee091b6aeac1113e09..638da99168860185f96753c917755b5c7a1b1a54 100644
--- a/examples/freeflowchannel/doc/_intro.md
+++ b/examples/freeflowchannel/doc/_intro.md
@@ -39,7 +39,7 @@ have a look at the DuMu<sup>x</sup> [documentation](https://dumux.org/docs/doxyg
 ## Problem set-up
 This example considers stationary flow of a fluid between two parallel solid plates in two dimensions.
 Flow is enforced from left to right by prescribing an inflow velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$
-on the left boundary, while a fixed pressure of $`p = 1.1 \text{bar}`$ and a zero velocity gradient
+on the left boundary, while a fixed pressure of $`p = 1.1 \, \text{bar}`$ and a zero velocity gradient
 in x-direction are prescribed on the right boundary. On the top and bottom boundaries, no-slip
 conditions are applied, which cause a parabolic velocity profile to develop along the channel.
 Take a look at Figure 1 for an illustration of the domain and the boundary conditions.
diff --git a/examples/freeflowchannel/main.cc b/examples/freeflowchannel/main.cc
index 2d2183d48a4b9b908cb5fe81904d06c7f722fcc5..b211d47995f3b28d78ce75bd139b4fea09cd7570 100644
--- a/examples/freeflowchannel/main.cc
+++ b/examples/freeflowchannel/main.cc
@@ -17,7 +17,7 @@
 #include <iostream>
 // [[/exclude]]
 
-// These are DUNE helper classes related to parallel computations and file I/O
+// These are DUNE helper classes related to parallel computations and file I/O.
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/grid/io/file/dgfparser/dgfexception.hh>
 
@@ -74,7 +74,7 @@ int main(int argc, char** argv) try
     Parameters::init(argc, argv);
     // [[/codeblock]]
 
-    // We define a convenience alias for the type tag of the problem. The type
+    // We define a convenience alias for the type tag of the problem (`Properties::TTag::ChannelExample`). The type
     // tag contains all the properties that are needed to define the model and the problem
     // setup. Throughout the main file, we will obtain types defined for this type tag
     // using the property system, i.e. with `GetPropType`.
@@ -94,16 +94,16 @@ int main(int argc, char** argv) try
 
     // #### Step 2: Setting up and solving the problem
     // First, a finite volume grid geometry is constructed from the grid that was created above.
-    // This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
+    // This builds the sub-control volumes (`scv`) and sub-control volume faces (`scvf`) for each element
     // of the grid partition.
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
     auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
 
-    // We now instantiate the problem, in which we define the boundary and initial conditions.
+    // We now instantiate the `problem`, in which we define the boundary and initial conditions.
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     auto problem = std::make_shared<Problem>(gridGeometry);
 
-    // We set a solution vector which consist of two parts: one part (indexed by `cellCenterIdx`)
+    // We set a solution vector `x` which consist of two parts: one part (indexed by `cellCenterIdx`)
     // is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
     // (indexed by `faceIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
     // We initialize the solution vector by what was defined as the initial solution of the the problem.
@@ -114,11 +114,11 @@ int main(int argc, char** argv) try
     problem->applyInitialSolution(x);
 
     // The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
-    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using GridVariables =GetPropType<TypeTag, Properties::GridVariables>;
     auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
     gridVariables->init(x);
 
-    // We then initialize the predefined model-specific output vtk output.
+    // We then initialize the predefined model-specific output VTK output.
     using IOFields = GetPropType<TypeTag, Properties::IOFields>;
     StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
     IOFields::initOutputModule(vtkWriter); // Add model specific output fields
@@ -126,13 +126,13 @@ int main(int argc, char** argv) try
 
     // [[details]] calculation of surface fluxes
     // We set up two surfaces over which fluxes are calculated.
-    // We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain.
-    // The first surface (added by the first call of addSurface) shall be placed at the middle of the channel.
+    // We determine the extend $`[xMin,xMax] \times [yMin,yMax]`$ of the physical domain.
+    // The first surface (added by the first call of `addSurface`) shall be placed at the middle of the channel.
     // If we have an odd number of cells in x-direction, there would not be any cell faces
     // at the position of the surface (which is required for the flux calculation).
     // In this case, we add half a cell-width to the x-position in order to make sure that
     // the cell faces lie on the surface. This assumes a regular cartesian grid.
-    // The second surface (second call of addSurface) is placed at the outlet of the channel.
+    // The second surface (second call of `addSurface`) is placed at the outlet of the channel.
     FluxOverSurface<GridVariables,
                     SolutionVector,
                     GetPropType<TypeTag, Properties::ModelTraits>,
@@ -166,7 +166,7 @@ int main(int argc, char** argv) try
     flux.addSurface("outlet", p0outlet, p1outlet);
     // [[/details]]
 
-    // We create and initialize the assembler for the stationary problem.
+    // We create and initialize the `assembler` for the stationary problem.
     // This is where the Jacobian matrix for the Newton solver is assembled.
     using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
diff --git a/examples/freeflowchannel/problem.hh b/examples/freeflowchannel/problem.hh
index 59b5cf27cbcb2f0767730c89d2ffb6f75531cc13..8d2d873fcacfac70ea7f7fb10169e4e706d5e286 100644
--- a/examples/freeflowchannel/problem.hh
+++ b/examples/freeflowchannel/problem.hh
@@ -24,7 +24,7 @@
 #include <dumux/freeflow/navierstokes/boundarytypes.hh>
 
 // ### The problem class
-// We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+// We enter the problem class `ChannelExampleProblem` where all necessary boundary conditions and initial conditions are set for our simulation.
 // As we are solving a problem related to free flow, we inherit from the base class `NavierStokesStaggeredProblem`.
 // [[codeblock]]
 namespace Dumux {
diff --git a/examples/freeflowchannel/properties.hh b/examples/freeflowchannel/properties.hh
index 08a022f0dd948b82c421b5fc8a443e77fbacc0b3..64be108beb1cc25afc1a0d5c525e53fafac2fb7f 100644
--- a/examples/freeflowchannel/properties.hh
+++ b/examples/freeflowchannel/properties.hh
@@ -10,7 +10,7 @@
 
 // ## Compile-time settings (`properties.hh`)
 //
-// In this file, the type tag used for this simulation is defined,
+// In this file, the type tag used for this simulation (`TTag::ChannelExample`) is defined,
 // for which we then specialize properties (compile time options) to the needs of the desired setup.
 //
 // [[content]]
@@ -19,7 +19,7 @@
 // [[details]] includes
 //
 // The `NavierStokes` type tag specializes most of the properties required for Navier-
-// Stokes single-phase flow simulations in DuMuX. We will use this in the following to inherit the
+// Stokes single-phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
 // respective properties and subsequently specialize those properties for our
 // type tag, which we want to modify or for which no meaningful default can be set.
 #include <dumux/freeflow/navierstokes/model.hh>
@@ -51,7 +51,7 @@
 // `Problem.EnableInertiaTerms = false`. Have a look at the input file `params.input`
 // to see how this is done in this example.
 // [[codeblock]]
-// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
+// We enter the namespace `Dumux::Properties` in order to import the entire `Dumux` namespace for general use:
 namespace Dumux::Properties {
 
 // declaration of the `ChannelExample` type tag for the single-phase flow problem
@@ -69,11 +69,11 @@ struct ChannelExample { using InheritsFrom = std::tuple<NavierStokes, StaggeredF
 template<class TypeTag>
 struct Grid<TypeTag, TTag::ChannelExample> { using type = Dune::YaspGrid<2>; };
 
-// This sets our problem class (see problem.hh) containing initial and boundary conditions.
+// This sets our problem type (see `problem.hh`) containing the initial and boundary conditions.
 template<class TypeTag>
 struct Problem<TypeTag, TTag::ChannelExample> { using type = Dumux::ChannelExampleProblem<TypeTag> ; };
 
-// This sets the fluid system to be used. Here, we use a liquid with constant properties as fluid phase.
+// This sets the fluid system type to be used. Here, we use a liquid with constant properties as the fluid phase.
 template<class TypeTag>
 struct FluidSystem<TypeTag, TTag::ChannelExample>
 {
@@ -86,7 +86,7 @@ struct FluidSystem<TypeTag, TTag::ChannelExample>
 // throughout the simulation.
 // [[details]] caching properties
 //
-// In Dumux, one has the option to activate/deactivate the grid-wide caching of
+// In DuMu<sup>x</sup>, one has the option to activate/deactivate the grid-wide caching of
 // geometries and variables. If active, the CPU time can be significantly reduced
 // as less dynamic memory allocation procedures are necessary. Per default, grid-wide
 // caching is disabled to ensure minimal memory requirements, however, in this example we
diff --git a/examples/liddrivencavity/doc/postprocessing.md b/examples/liddrivencavity/doc/postprocessing.md
index 600b0d57ec667c5085dac6e9ad9353c48cb88672..9c88c88c35ec43f97697f6036611e9c3d89c8485 100644
--- a/examples/liddrivencavity/doc/postprocessing.md
+++ b/examples/liddrivencavity/doc/postprocessing.md
@@ -5,7 +5,7 @@
 |---|---:|
 
 # Part 2: Post processing
-In this part we first visualize our simulation result and then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
+In this part we first visualize our simulation result. Then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
 
 ## Step 1. Visualize results with Paraview
 
@@ -27,7 +27,7 @@ For Re = 1 and Re = 1000, the result should look like this:
 ## Step 2. Compare our data with reference
 
 The verification and validation are essential to guarantee the accuracy and credibility of the numerical models.
-The velocity components for the velocity components at x = 0.5m and y = 0.5m are obtained as we run the test cases. 
+The velocity components for the velocity components at x = 0.5~m and y = 0.5~m are obtained as we run the test cases. 
 
 We compare our results with the reference data reconstructed from [Ghia et al.](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
 For convenience, we placed the reference data in the folder named `reference_data`. For instance, the files `ghia_x.csv` and `ghia_y.csv` represent the reference vertical velocity component $`v_x`$ at x = 0.5 m and $`v_y`$ at y = 0.5 m for the scenario Re = 1000, respectively.
diff --git a/examples/liddrivencavity/doc/postprocessing_text.md b/examples/liddrivencavity/doc/postprocessing_text.md
index 02cc22c670fd90865dd64eeab15b85c743bd356c..b71764a7595fb098e1b8dcd5feed778d02e9fa46 100644
--- a/examples/liddrivencavity/doc/postprocessing_text.md
+++ b/examples/liddrivencavity/doc/postprocessing_text.md
@@ -1,5 +1,5 @@
 # Part 2: Post processing
-In this part we first visualize our simulation result and then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
+In this part we first visualize our simulation result. Then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
 
 ## Step 1. Visualize results with Paraview
 
@@ -21,7 +21,7 @@ For Re = 1 and Re = 1000, the result should look like this:
 ## Step 2. Compare our data with reference
 
 The verification and validation are essential to guarantee the accuracy and credibility of the numerical models.
-The velocity components for the velocity components at x = 0.5m and y = 0.5m are obtained as we run the test cases. 
+The velocity components for the velocity components at x = 0.5~m and y = 0.5~m are obtained as we run the test cases. 
 
 We compare our results with the reference data reconstructed from [Ghia et al.](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
 For convenience, we placed the reference data in the folder named `reference_data`. For instance, the files `ghia_x.csv` and `ghia_y.csv` represent the reference vertical velocity component $`v_x`$ at x = 0.5 m and $`v_y`$ at y = 0.5 m for the scenario Re = 1000, respectively.
diff --git a/examples/liddrivencavity/doc/problem.md b/examples/liddrivencavity/doc/problem.md
index a41a2b43c38908a951f54dc7682ea39d9ca542c3..86380bd571aacbe8994c0a3567389e36f1a86085 100644
--- a/examples/liddrivencavity/doc/problem.md
+++ b/examples/liddrivencavity/doc/problem.md
@@ -179,7 +179,7 @@ Include the `NavierStokesBoundaryTypes` class which specifies the boundary types
 
 ### The problem class
 As we are solving a problem related to free flow, we create a new class called `LidDrivenCavityExampleProblem`
-and let it inherit from a base class for the momentum and mass subproblems (selected in properties.hh).
+and let it inherit from a base class for the momentum and mass subproblems (selected in `properties.hh`).
 
 ```cpp
 namespace Dumux {
@@ -334,7 +334,7 @@ The following function defines the initial conditions.
     }
 ```
 
-the data members of the problem class
+Finally, the (private) data members of the problem class.
 
 ```cpp
 private:
@@ -604,7 +604,7 @@ the non-linear solver
 
 ##### The time loop
 In each time step, we solve the non-linear system of equations, write
-the current solution into .vtk files and prepare for the next time step.
+the current solution into VTK files and prepare for the next time step.
 
 ```cpp
     timeLoop->start(); do
@@ -632,7 +632,7 @@ the current solution into .vtk files and prepare for the next time step.
     } while (!timeLoop->finished());
 ```
 
-We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file
+We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file.
 
 ```cpp
     writeSteadyVelocityAndCoordinates(*momentumProblem, x[momentumIdx]);
diff --git a/examples/liddrivencavity/main.cc b/examples/liddrivencavity/main.cc
index de533c79c4a857eb06c9a2fae0a05cda1d4bb9c5..8a9275c16702e26fbf247044551a3654926a1e09 100644
--- a/examples/liddrivencavity/main.cc
+++ b/examples/liddrivencavity/main.cc
@@ -206,7 +206,7 @@ int main(int argc, char** argv)
 
     // ##### The time loop
     // In each time step, we solve the non-linear system of equations, write
-    // the current solution into .vtk files and prepare for the next time step.
+    // the current solution into VTK files and prepare for the next time step.
     // [[codeblock]]
     timeLoop->start(); do
     {
@@ -233,7 +233,7 @@ int main(int argc, char** argv)
     } while (!timeLoop->finished());
     // [[/codeblock]]
 
-    // We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file
+    // We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file.
     writeSteadyVelocityAndCoordinates(*momentumProblem, x[momentumIdx]);
 
     // The following piece of code prints a final status report of the time loop
diff --git a/examples/liddrivencavity/problem.hh b/examples/liddrivencavity/problem.hh
index 8798efac2d402166c31783c220095ddd19d0ffb9..b90d76fb1d55fe468654cdaffcb50e60ef09aed1 100644
--- a/examples/liddrivencavity/problem.hh
+++ b/examples/liddrivencavity/problem.hh
@@ -25,7 +25,7 @@
 
 // ### The problem class
 // As we are solving a problem related to free flow, we create a new class called `LidDrivenCavityExampleProblem`
-// and let it inherit from a base class for the momentum and mass subproblems (selected in properties.hh).
+// and let it inherit from a base class for the momentum and mass subproblems (selected in `properties.hh`).
 // [[codeblock]]
 namespace Dumux {
 template <class TypeTag, class BaseProblem>
@@ -172,7 +172,7 @@ public:
         return values;
     }
     // [[/codeblock]]
-    // the data members of the problem class
+    // Finally, the (private) data members of the problem class.
     // [[codeblock]]
 private:
     static constexpr Scalar eps_ = 1e-6;
diff --git a/examples/porenetwork_upscaling/README.md b/examples/porenetwork_upscaling/README.md
index 50076569d5e19463beaebbf7653fa137b9777387..b59f1a5d4a113bb4e8d6128b8097b76ee12084b3 100644
--- a/examples/porenetwork_upscaling/README.md
+++ b/examples/porenetwork_upscaling/README.md
@@ -6,12 +6,10 @@ __In this example, you will learn how to__
 
 * simulate creeping/non-creeping flow on a pore network by applying a pressure gradient in a given direction
 * perform an upscaling in order to determine the flow properties of the porous medium such as:
-the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using the creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using the non-creeping flow simulation.
-
+the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using a creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using a non-creeping flow simulation.
 
 __Result__.
-As a result of the creeping flow simulation of this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction:
-
+As a result of the creeping flow simulation in this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction:
 ```
 X-direction:
 -- Darcy (intrinsic) permeability = 3.326e-12 m^2
@@ -28,7 +26,7 @@ Figure 1 shows the pressure distribution within the pore network for the case of
     </center>
 </figure>
 
-The non-creeping flow simulation additionally gives you the Forchheimer permeability and coefficient for each spatial direction as in the example below for the x-direction:
+The non-creeping flow simulation additionally gives you the Forchheimer permeabilities and coefficients for each spatial direction as in the example below for the x-direction:
 
 ```
 X-direction:
diff --git a/examples/porenetwork_upscaling/doc/_intro.md b/examples/porenetwork_upscaling/doc/_intro.md
index 376f3c0363c1880adc08cea0a8b7a93bf40623e2..d83502f45e95d75ae0ce4eb996182843e76956ba 100644
--- a/examples/porenetwork_upscaling/doc/_intro.md
+++ b/examples/porenetwork_upscaling/doc/_intro.md
@@ -4,12 +4,10 @@ __In this example, you will learn how to__
 
 * simulate creeping/non-creeping flow on a pore network by applying a pressure gradient in a given direction
 * perform an upscaling in order to determine the flow properties of the porous medium such as:
-the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using the creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using the non-creeping flow simulation.
-
+the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using a creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using a non-creeping flow simulation.
 
 __Result__.
-As a result of the creeping flow simulation of this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction:
-
+As a result of the creeping flow simulation in this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction:
 ```
 X-direction:
 -- Darcy (intrinsic) permeability = 3.326e-12 m^2
@@ -26,7 +24,7 @@ Figure 1 shows the pressure distribution within the pore network for the case of
     </center>
 </figure>
 
-The non-creeping flow simulation additionally gives you the Forchheimer permeability and coefficient for each spatial direction as in the example below for the x-direction:
+The non-creeping flow simulation additionally gives you the Forchheimer permeabilities and coefficients for each spatial direction as in the example below for the x-direction:
 
 ```
 X-direction:
diff --git a/examples/porenetwork_upscaling/doc/main.md b/examples/porenetwork_upscaling/doc/main.md
index 5d646066be7573fdda261ec0081a0986eee98053..8d3c35b447d531d8a4c5b078c723c60005150dda 100644
--- a/examples/porenetwork_upscaling/doc/main.md
+++ b/examples/porenetwork_upscaling/doc/main.md
@@ -108,9 +108,9 @@ void runExample()
 
 ### Instantiate the solver
 We use the `NewtonSolver` class, which is instantiated on the basis
-of an assembler and a linear solver. When the `solve` function of the
-`NewtonSolver` is called, it uses the assembler and linear
-solver classes to assemble and solve the non-linear system.
+of an assembler and a linear solver. When the `solve()` function of the
+Newton solver instance is called, it uses the assembler and
+linear solver to assemble and solve the linear systems in each Newton step.
 
 ```cpp
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
diff --git a/examples/porenetwork_upscaling/doc/problem.md b/examples/porenetwork_upscaling/doc/problem.md
index f6d1871a922662bd686e1d7e1d5a8d510459cfbf..7d8144665d7245a1fd6cfc0ec63c7c5be94b82aa 100644
--- a/examples/porenetwork_upscaling/doc/problem.md
+++ b/examples/porenetwork_upscaling/doc/problem.md
@@ -76,8 +76,8 @@ The classes that define the problem and parameters used in this simulation
 </details>
 
 ### `TypeTag` definition
-Two `TypeTag` for our simulation are defined, one for creeping flow and another for non-creeping flow,
-which inherit properties from the single-phase pore network model. The non-creeping flow inherits
+Two type tags for our simulation are defined, one for creeping flow (`PNMUpscalingCreepingFlow`) and another for non-creeping flow (`PNMUpscalingNonCreepingFlow`),
+which inherit properties from the single-phase pore network model (`ONMOneP`). The non-creeping flow inherits
 all properties from the creeping flow simulation but overrides the `AdvectionType` property.
 
 ```cpp
diff --git a/examples/porenetwork_upscaling/doc/upscalinghelper.md b/examples/porenetwork_upscaling/doc/upscalinghelper.md
index 377b6493311fc12c90da43a5ec04cece9e8e9824..f6b7a07e7980b016e49f3b63c1b732360147516a 100644
--- a/examples/porenetwork_upscaling/doc/upscalinghelper.md
+++ b/examples/porenetwork_upscaling/doc/upscalinghelper.md
@@ -4,7 +4,7 @@
 | [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 2](main.md) |
 |---|---:|
 
-The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in this direction. Firstly, it evaluates the apparent velocity as:
+The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in a chosen direction. Firstly, it evaluates the apparent velocity as:
 
 ```math
      v_{\mathrm{Apparent},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}}
@@ -32,7 +32,7 @@ To calculate upscaled properties, we rearrange Forchehimer's equation as:
 ```math
  \frac{\nabla p_i v_{\mathrm{Apparent},i}}{\mu} = \frac{1}{K_f} + \frac{\varrho v_{\mathrm{Apparent},i}}{\mu} \beta .
 ```
-Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in params.input. it is recommended to use more than 10 pressure sample points which can be set in the input file. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately.
+Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in `params.input`. It is recommended to use more than 10 pressure sample points which can be set in the input file using `Problem.NumberOfPressureGradients`. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately.
 
 The code documentation is structured as follows:
 
diff --git a/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md b/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md
index 76c143622b1edb753a7bd4fd153ffafb0ace4b63..66f5ac5e910aa7295b949c5f1fd3d1a55cd5d540 100644
--- a/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md
+++ b/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md
@@ -1,4 +1,4 @@
-The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in this direction. Firstly, it evaluates the apparent velocity as:
+The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in a chosen direction. Firstly, it evaluates the apparent velocity as:
 
 ```math
      v_{\mathrm{Apparent},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}}
@@ -26,6 +26,6 @@ To calculate upscaled properties, we rearrange Forchehimer's equation as:
 ```math
  \frac{\nabla p_i v_{\mathrm{Apparent},i}}{\mu} = \frac{1}{K_f} + \frac{\varrho v_{\mathrm{Apparent},i}}{\mu} \beta .
 ```
-Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in params.input. it is recommended to use more than 10 pressure sample points which can be set in the input file. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately.
+Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in `params.input`. It is recommended to use more than 10 pressure sample points which can be set in the input file using `Problem.NumberOfPressureGradients`. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately.
 
 The code documentation is structured as follows:
diff --git a/examples/porenetwork_upscaling/main.cc b/examples/porenetwork_upscaling/main.cc
index e58ea99b662b5af9caa1e2f62eb6b9a07649293f..1f16dc40900cdb8d1c96f8a2960bd3dada0bbea5 100644
--- a/examples/porenetwork_upscaling/main.cc
+++ b/examples/porenetwork_upscaling/main.cc
@@ -86,9 +86,9 @@ void runExample()
 
     // ### Instantiate the solver
     // We use the `NewtonSolver` class, which is instantiated on the basis
-    // of an assembler and a linear solver. When the `solve` function of the
-    // `NewtonSolver` is called, it uses the assembler and linear
-    // solver classes to assemble and solve the non-linear system.
+    // of an assembler and a linear solver. When the `solve()` function of the
+    // Newton solver instance is called, it uses the assembler and
+    // linear solver to assemble and solve the linear systems in each Newton step.
     // [[codeblock]]
     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
     auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
diff --git a/examples/porenetwork_upscaling/properties.hh b/examples/porenetwork_upscaling/properties.hh
index 46111ffe8c6751fdb2e32a6a921d23919b3f190c..f33a880cc7a70e3b8f94582e52bc5dc078287b3f 100644
--- a/examples/porenetwork_upscaling/properties.hh
+++ b/examples/porenetwork_upscaling/properties.hh
@@ -47,8 +47,8 @@
 // [[/details]]
 //
 // ### `TypeTag` definition
-// Two `TypeTag` for our simulation are defined, one for creeping flow and another for non-creeping flow,
-// which inherit properties from the single-phase pore network model. The non-creeping flow inherits
+// Two type tags for our simulation are defined, one for creeping flow (`PNMUpscalingCreepingFlow`) and another for non-creeping flow (`PNMUpscalingNonCreepingFlow`),
+// which inherit properties from the single-phase pore network model (`ONMOneP`). The non-creeping flow inherits
 // all properties from the creeping flow simulation but overrides the `AdvectionType` property.
 namespace Dumux::Properties {
 namespace TTag {
diff --git a/examples/shallowwaterfriction/README.md b/examples/shallowwaterfriction/README.md
index 7ed5de42250b74650c28ce296fd132242dd73126..f2cc22d2bc15213828a0aea203e7d98ebe813668 100644
--- a/examples/shallowwaterfriction/README.md
+++ b/examples/shallowwaterfriction/README.md
@@ -25,8 +25,8 @@ __Table of contents__. This description is structured as follows:
 ## Problem set-up
 ### Model domain
 The model domain is given by a rough channel with a slope of 0.001.
-The domain is 500 meters long and 5 meters wide.
-The bottom altitude is 10 m at the inflow and hence 9.5 m at the outflow.
+The domain is $`500 \, \mathrm{m}`$ long and $`5 mathrm{m}`$ wide.
+The bottom altitude is $`10 mathrm{m}`$ at the inflow and hence $`9.5 mathrm{m}`$ at the outflow.
 Bottom friction is considered by applying
 [Manning's law](#mannings-law) ($`n`$ = 0.025).
 
@@ -35,15 +35,15 @@ At the lateral sides a no-flow boundary condition is applied. Also, no friction
 considered there and therefore a no slip boundary
 condition is applied. These are the default boundary conditions for the shallow
 water model. At the left border, a discharge boundary condition
-is applied as inflow boundary condition with $`q = -1.0 m^2 s^{-1}`$.
+is applied as inflow boundary condition with $`q = -1.0 mathrm{m}^2 \mathrm{s}^{-1}`$.
 At the right border, a fixed water depth boundary condition
 is applied for the outflow. Normal flow is assumed, therefore the water
 depth at the right border is calculated using the equation
 of [Gauckler, Manning and Strickler](#analytical-solution).
 
 ### Initial conditions
-The initial water depth is set to 1 m, which is slightly higher than the normal flow
-water depth (0.87 m). Therefore, we expect a decreasing
+The initial water depth is set to 1~m, which is slightly higher than the normal flow
+water depth (0.87~m). Therefore, we expect a decreasing
 water level during the simulation until the normal flow condition is reached in
 the entire model domain. The initial velocity is set to zero.
 
diff --git a/examples/shallowwaterfriction/doc/_intro.md b/examples/shallowwaterfriction/doc/_intro.md
index d5328920d47a4c9ad1989988d5d0daa1727b92bf..bd1e7d2e8b4f8db4a70e93befbe4af6b2b62e8b2 100644
--- a/examples/shallowwaterfriction/doc/_intro.md
+++ b/examples/shallowwaterfriction/doc/_intro.md
@@ -23,8 +23,8 @@ __Table of contents__. This description is structured as follows:
 ## Problem set-up
 ### Model domain
 The model domain is given by a rough channel with a slope of 0.001.
-The domain is 500 meters long and 5 meters wide.
-The bottom altitude is 10 m at the inflow and hence 9.5 m at the outflow.
+The domain is $`500 \, \mathrm{m}`$ long and $`5 mathrm{m}`$ wide.
+The bottom altitude is $`10 mathrm{m}`$ at the inflow and hence $`9.5 mathrm{m}`$ at the outflow.
 Bottom friction is considered by applying
 [Manning's law](#mannings-law) ($`n`$ = 0.025).
 
@@ -33,15 +33,15 @@ At the lateral sides a no-flow boundary condition is applied. Also, no friction
 considered there and therefore a no slip boundary
 condition is applied. These are the default boundary conditions for the shallow
 water model. At the left border, a discharge boundary condition
-is applied as inflow boundary condition with $`q = -1.0 m^2 s^{-1}`$.
+is applied as inflow boundary condition with $`q = -1.0 mathrm{m}^2 \mathrm{s}^{-1}`$.
 At the right border, a fixed water depth boundary condition
 is applied for the outflow. Normal flow is assumed, therefore the water
 depth at the right border is calculated using the equation
 of [Gauckler, Manning and Strickler](#analytical-solution).
 
 ### Initial conditions
-The initial water depth is set to 1 m, which is slightly higher than the normal flow
-water depth (0.87 m). Therefore, we expect a decreasing
+The initial water depth is set to 1~m, which is slightly higher than the normal flow
+water depth (0.87~m). Therefore, we expect a decreasing
 water level during the simulation until the normal flow condition is reached in
 the entire model domain. The initial velocity is set to zero.