doc: remove unicode characters, other fixes

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
 <img src="real-sandbox-1.jpg" style="width:45.0%" />
 
- Examples of deformation patterns of &#x201C;sand box&#x201D; 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
 <img src="real-sandbox-2.jpg" style="width:45.0%" />
 
- Examples of deformation patterns of &#x201C;sand box&#x201D; 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

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
 <img src="crameri_2_convergence.svg" style="width:100.0%" />
 
-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 <span class="math inline">&#x2004;&#x223C;&#x2004;0.1</span>, 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$

benchmarks/gravity_mantle/doc/gravity_mantle.md

@@ -91,11 +91,11 @@ information.
 ```{figure-md} fig:grav_mantle4
 <img src="grav.*" style="width:48.0%" />
 
- Mantle gravity: gravitational acceleration |g| computed at radius 6621&#xA0;km.
+ Mantle gravity: gravitational acceleration |g| computed at radius $6621~\text{ km}$.
 ```
 
 ```{figure-md} fig:grav_mantle5
 <img src="pot.*" style="width:48.0%" />
 
- Mantle gravity: gravitational potential computed at radius 6621&#xA0;km.
+ Mantle gravity: gravitational potential computed at radius $6621~\text{ km}$.
 ```

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
-<img src="inclusion-solution.*" alt="The viscosity field when interpolated onto the mesh (internally, the &#x201C;exact&#x201D; viscosity field &#x2013; large inside a circle, small outside &#x2013; is used), and overlaid to it some velocity vectors." />
+<img src="inclusion-solution.*" />
 
-The viscosity field when interpolated onto the mesh (internally, the &#x201C;exact&#x201D; viscosity field &#x2013; large inside a circle, small outside &#x2013; 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
-<img src="inclusion-solution-pressure.*"  alt="The pressure with its oscillations along the interface. The oscillations become more localized as the mesh is refined." />
+<img src="inclusion-solution-pressure.*" />
 
 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).

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

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&nbsp;m<sup>&minus;3</sup>, $\alpha = 2.5\cdot10^{-5}$
-&nbsp;K, $k = 4$ W&nbsp;m<sup>&minus;1</sup>&nbsp;K<sup>&minus;1</sup>, $C_p = 1000$
-J&nbsp;kg<sup>&minus;1</sup>&nbsp;K<sup>&minus;1</sup> and its internal heating rate is
-$Q{_{int}} = 1\cdot10^{-12}$ W&nbsp;kg<sup>&minus;1</sup>. 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
 <img src="temps.png" style="width:90.0%" />
 
- Bunge et al.&#xA0;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.
 ```

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
 <img src="box_3d_mesh.*" style="width:100.0%" />
 
-Convection in a 3d box: Meshes and large-scale velocity field for the third, fourth and sixth of the snapshots shown in Fig.&#xA0;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

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

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

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
 <img src="morency_doin_2004_fig1.svg" />
 
-Approximate reproduction of figure 1 from <span class="citation" data-cites="MD04">{cite:t}`morency:doin:2004`</span> using the &#x2018;morency doin&#x2019; 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 <span class="citation" data-cites="MD04">{cite:t}`morency:doin:2004`</span> 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.
 ```

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
 <img src="init_diff_rms_vel_comparison.png" style="width:50.0%" />
 
- 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&#xA0;[fig:vk-6] and&#xA0;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