diff --git a/benchmarks/buiter_et_al_2016_jsg/doc/buiter_et_al_2016_jsg.md b/benchmarks/buiter_et_al_2016_jsg/doc/buiter_et_al_2016_jsg.md
index 5f7ca230bc3..deb2174f08b 100644
--- a/benchmarks/buiter_et_al_2016_jsg/doc/buiter_et_al_2016_jsg.md
+++ b/benchmarks/buiter_et_al_2016_jsg/doc/buiter_et_al_2016_jsg.md
@@ -17,13 +17,13 @@ the patterns one can then observe in these do-it-yourself models are shown in
```{figure-md} fig:sandbox-images1
- Examples of deformation patterns of “sand box” experiments in which alternating layers of differently-colored sand undergo deformation. Pictures courtesy of the lab of Dennis Harry at Colorado State University.
+ Examples of deformation patterns of sand box experiments in which alternating layers of differently-colored sand undergo deformation. Pictures courtesy of the lab of Dennis Harry at Colorado State University.
```
```{figure-md} fig:sandbox-images2
- Examples of deformation patterns of “sand box” experiments in which alternating layers of differently-colored sand undergo deformation. Pictures courtesy of the lab of Dennis Harry at Colorado State University.
+ Examples of deformation patterns of sand box experiments in which alternating layers of differently-colored sand undergo deformation. Pictures courtesy of the lab of Dennis Harry at Colorado State University.
```
{cite}`buiter:etal:2016` organized new comparison experiments
diff --git a/benchmarks/crameri_et_al/doc/crameri_et_al.md b/benchmarks/crameri_et_al/doc/crameri_et_al.md
index cfc8da280c6..e4553f716bb 100644
--- a/benchmarks/crameri_et_al/doc/crameri_et_al.md
+++ b/benchmarks/crameri_et_al/doc/crameri_et_al.md
@@ -102,7 +102,7 @@ CFL number and mesh resolution. The results are shown in {numref}`fig:crameri-be
```{figure-md} fig:crameri-benchmark-convergence
-Convergence for case two. Left: Logarithm of the error with decreasing CFL number. As the CFL number decreases, the error gets smaller. However, once it reaches a value of ∼ 0.1, there stops being much improvement in accuracy. Right: Logarithm of the error with increasing maximum mesh resolution. As the resolution increases, so does the accuracy.
+Convergence for case two. Left: Logarithm of the error with decreasing CFL number. As the CFL number decreases, the error gets smaller. However, once it reaches a value of $\sim 0.1$, there stops being much improvement in accuracy. Right: Logarithm of the error with increasing maximum mesh resolution. As the resolution increases, so does the accuracy.
```
We find that at 3 Ma converges to a maximum topography of $\sim 396$
diff --git a/benchmarks/gravity_mantle/doc/gravity_mantle.md b/benchmarks/gravity_mantle/doc/gravity_mantle.md
index 6114863f529..3a0d8c40add 100644
--- a/benchmarks/gravity_mantle/doc/gravity_mantle.md
+++ b/benchmarks/gravity_mantle/doc/gravity_mantle.md
@@ -91,11 +91,11 @@ information.
```{figure-md} fig:grav_mantle4
- Mantle gravity: gravitational acceleration |g| computed at radius 6621 km.
+ Mantle gravity: gravitational acceleration |g| computed at radius $6621~\text{ km}$.
```
```{figure-md} fig:grav_mantle5
- Mantle gravity: gravitational potential computed at radius 6621 km.
+ Mantle gravity: gravitational potential computed at radius $6621~\text{ km}$.
```
diff --git a/benchmarks/inclusion/doc/inclusion.md b/benchmarks/inclusion/doc/inclusion.md
index 00d99378f0b..aa27e6e4f83 100644
--- a/benchmarks/inclusion/doc/inclusion.md
+++ b/benchmarks/inclusion/doc/inclusion.md
@@ -13,15 +13,17 @@ numerical solution. This can be seen in the visualizations shown in
analytic solution against which we compare is given in {cite:t}`schmid:podladchikov:2003`. An extensive discussion of convergence properties is given in {cite:t}`kronbichler:etal:2012`.
```{figure-md} fig:inclusion1
-
+
-The viscosity field when interpolated onto the mesh (internally, the “exact” viscosity field – large inside a circle, small outside – is used), and overlaid to it some velocity vectors.
+The viscosity field when interpolated onto the mesh (internally, the exact viscosity field is used): large inside a circle and small outside it. Overlaid to it a vector field showing the velocity.
```
+
```{figure-md} fig:inclusion2
-
+
The pressure with its oscillations along the interface. The oscillations become more localized as the mesh is refined.
```
+
The benchmark can be run using the parameter files in
[benchmarks/inclusion/](https://github.com/geodynamics/aspect/blob/main/benchmarks/inclusion). The material model, boundary condition, and
postprocessor are defined in [benchmarks/inclusion/inclusion.cc](https://github.com/geodynamics/aspect/blob/main/benchmarks/inclusion/inclusion.cc).
diff --git a/benchmarks/rayleigh_taylor_instability/doc/rayleigh_taylor_instability.md b/benchmarks/rayleigh_taylor_instability/doc/rayleigh_taylor_instability.md
index 9a570a90d3a..8afd0970dc1 100644
--- a/benchmarks/rayleigh_taylor_instability/doc/rayleigh_taylor_instability.md
+++ b/benchmarks/rayleigh_taylor_instability/doc/rayleigh_taylor_instability.md
@@ -40,10 +40,10 @@ j_{22} &= \frac{\eta_1 2 \phi_1^2 \phi_2}{\eta_2(\cosh 2\phi_1 -1-2\phi_1^2)} -
\phi_2&=\frac{2\pi h_2}{\lambda}
```
We set
-$L_x=L_y=\text{ 512 km}$, $h_1=h_2=\text{ 256 km}$,
-$|\boldsymbol{g}|=\text{10 m/s^2}$, $\Delta=\text{3 km}$,
-$\rho_1=\text{3300 kg/m^3}$, $\rho_2=\text{3000 kg/m^3}$,
-$\eta_1=\text{1e21 Pa.s}$. $\eta_2$ is varied between $10^{20}$ and $10^{23}$
+$L_x=L_y=512 \text{ km}$, $h_1=h_2=256 \text{ km}$,
+$|\boldsymbol{g}|=10\text{m}/\text{s}^2$, $\Delta=3\text{ km}$,
+$\rho_1=3300\text{ kg}/\text{m}^3$, $\rho_2=3000\text{kg}/\text{m}^3$,
+$\eta_1=1e21\text{ Pa s}$. $\eta_2$ is varied between $10^{20}$ and $10^{23}$
and 3 values of $\lambda$ (64, 128, and 256km) are used. Adaptive mesh
refinement based on density is used to capture the interface between the two
fluids, as shown in {numref}`fig:RTi_grids_a` and {numref}`fig:RTi_grids_b`. This translates as follows in the input
diff --git a/cookbooks/bunge_et_al_mantle_convection/doc/bunge_et_al_mantle_convection.md b/cookbooks/bunge_et_al_mantle_convection/doc/bunge_et_al_mantle_convection.md
index 2f5bc67dc41..cc75ee2c489 100644
--- a/cookbooks/bunge_et_al_mantle_convection/doc/bunge_et_al_mantle_convection.md
+++ b/cookbooks/bunge_et_al_mantle_convection/doc/bunge_et_al_mantle_convection.md
@@ -21,10 +21,8 @@ $T{_{cmb}} = 3450$ K. The gravity vector is radial and its
magnitude is $g = 10 \text{m s}^{-2}$.
There is a single incompressible fluid in the domain, characterized by
-$\rho_0 = 4500$ kg m−3, $\alpha = 2.5\cdot10^{-5}$
- K, $k = 4$ W m−1 K−1, $C_p = 1000$
-J kg−1 K−1 and its internal heating rate is
-$Q{_{int}} = 1\cdot10^{-12}$ W kg−1. The
+$\rho_0 = 4500\text{ kg m}^{-3}$, $\alpha = 2.5\cdot10^{-5}\text{ K}$, $k = 4\text{ W m}^{-1}\text{K}^{-1}$, $C_p = 1000 \text{ J kg}^{-1}\text{K}^{-1}$ and its internal heating rate is
+$Q{_{int}} = 1\cdot10^{-12}\text{ W kg}^{-1}$. The
interface between the upper mantle (viscosity $\eta_{um}$) and
the lower mantle (viscosity $\eta_{lm}$) is fixed at 670 km
depth. As in the article we consider four time-independent radial viscosity
@@ -56,5 +54,5 @@ movie of how the temperature evolves over this time period at
```{figure-md} fig:bunge_et_al
- Bunge et al. benchmark. From left to right: temperature field at time $t=5\cdot 10^9$ years obtained with viscosity profiles a, b, c and d.
+ Bunge et al. benchmark. From left to right: temperature field at time $t=5\cdot 10^9$ years obtained with viscosity profiles a, b, c and d.
```
diff --git a/cookbooks/convection_box_3d/doc/convection_box_3d.md b/cookbooks/convection_box_3d/doc/convection_box_3d.md
index 82904111ea0..5089e8251ca 100644
--- a/cookbooks/convection_box_3d/doc/convection_box_3d.md
+++ b/cookbooks/convection_box_3d/doc/convection_box_3d.md
@@ -140,7 +140,7 @@ Convection in a 3d box: Temperature isocontours and some velocity vectors at the
```{figure-md} fig:box-3d-mesh
-Convection in a 3d box: Meshes and large-scale velocity field for the third, fourth and sixth of the snapshots shown in Fig. 6.
+Convection in a 3d box: Meshes and large-scale velocity field for the third, fourth and sixth of the snapshots shown in Fig.{numref}`fig:box-3d-solution`.
```
As before, we could analyze all sorts of data from the statistics file but we
diff --git a/cookbooks/global_melt/doc/global_melt.md b/cookbooks/global_melt/doc/global_melt.md
index bee1c0f7059..19e0c6704c4 100644
--- a/cookbooks/global_melt/doc/global_melt.md
+++ b/cookbooks/global_melt/doc/global_melt.md
@@ -17,7 +17,7 @@ input file has to be modified in order to add melt transport. A movie that
compares the evolution of the temperature field and the amount of melt present
in both models in higher resolution can be found [online](https://www.youtube.com/watch?v=Kwyp4Jvx6MU).
-The model setup is a 2D box with dimensions of $2900 \times 8700$ km,
+The model setup is a 2D box with dimensions of $2900 \times 8700$ km,
and it is heated from the bottom and cooled from the top. A full description
can be found in Section 4.7 "Influence of melt migration on a
global-scale convection model" in {cite:t}`dannberg:heister:2016`. In the
diff --git a/cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md b/cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md
index 09eaf3e2c8f..fece83c5412 100644
--- a/cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md
+++ b/cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md
@@ -21,7 +21,7 @@ previous cookbook {ref}`sec:cookbooks:global-melt`.
As the flow at mid-ocean ridges can be assumed to be roughly symmetric with
respect to the ridge axis in the center, we only model one half of the ridge
-in a 2d Cartesian box with dimensions of $105 \times 70$ km. Solid
+in a 2d Cartesian box with dimensions of $105 \times 70$ km. Solid
material is flowing in from the bottom with a prescribed temperature and
melting due to decompression as is rises. The model is cooled from the top so
that melt freezes again as it approaches this boundary. In addition, a fixed
diff --git a/cookbooks/morency_doin_2004/doc/morency_doin_2004.md b/cookbooks/morency_doin_2004/doc/morency_doin_2004.md
index 4c313c942e2..132f14fa8c1 100644
--- a/cookbooks/morency_doin_2004/doc/morency_doin_2004.md
+++ b/cookbooks/morency_doin_2004/doc/morency_doin_2004.md
@@ -55,5 +55,5 @@ The viscosity profile in Figure 1 of {cite:t}`morency:doin:2004` appears to be w
```{figure-md} fig:md-1
-Approximate reproduction of figure 1 from {cite:t}`morency:doin:2004` using the ‘morency doin’ material model with almost all default parameters. Note the low-viscosity Moho, enabled by the low activation energy of the crustal component.
+Approximate reproduction of figure 1 from {cite:t}`morency:doin:2004` using the "morency doin" material model with almost all default parameters. Note the low-viscosity Moho, enabled by the low activation energy of the crustal component.
```
diff --git a/cookbooks/van-keken-vof/doc/van-keken-vof.md b/cookbooks/van-keken-vof/doc/van-keken-vof.md
index 77fa909cccb..46c23f3a40b 100644
--- a/cookbooks/van-keken-vof/doc/van-keken-vof.md
+++ b/cookbooks/van-keken-vof/doc/van-keken-vof.md
@@ -176,7 +176,7 @@ the cosine function in the initial conditions, as shown below.
```{figure-md} fig:vof-vk-3
- Computations of the van Keken problem made with the VOF interface tracking algorithm showing the evolution of the RMS velocity as a function of time for small changes in the amplitude a of the cosine function in the initial condition at 7 levels of refinement. Compare to Figures [fig:vk-6] and 1.
+ Computations of the van Keken problem made with the VOF interface tracking algorithm showing the evolution of the RMS velocity as a function of time for small changes in the amplitude a of the cosine function in the initial condition at 7 levels of refinement. Compare to {numref}`fig:vk-6` and {numref}`fig:vk-9`.
```
In these computations we vary the value of $a$ from its usual value of