diff --git a/Project.toml b/Project.toml
index 89907ffe..4a796a36 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GeophysicalModelGenerator"
 uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 authors = ["Boris Kaus", "Marcel Thielmann"]
-version = "0.7.9"
+version = "0.7.10"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
diff --git a/README.md b/README.md
index a524aec6..a13bc2cb 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@
 
 Creating consistent 3D images of geophysical and geological datasets and turning that into an input model for geodynamic simulations is often challenging. The aim of this package is to help with this, by providing a number of routines to easily import data and create a consistent 3D visualisation from it in the VTK-toolkit format, which can for example be viewed with [Paraview](https://www.paraview.org). In addition, we provide a range of tools that helps to generate input models to perform geodynamic simulations and import the results of such simulations back into julia.
 
-A short summary of the package and its features are given below. For a detailed description of the package and to learn how to use it, have a look at the [documentation](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev/). 
+A short summary of the package and its features are given below. For a detailed description of the package and to learn how to use it, have a look at the [documentation](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev/).
 
 ![README_img](./docs/src/assets/img/Readme_pic.png)
 ### Contents
@@ -43,17 +43,17 @@ The best way to learn how to use this is to install the package (see below) and
 ## Installation
 First, you need to install julia on your machine. We recommend to use the binaries from [https://julialang.org](https://julialang.org).
 Next, start julia and switch to the julia package manager using `]`, after which you can add the package.
-```julia
+```julia-repl
 julia> ]
-(@v1.10) pkg> add GeophysicalModelGenerator
+(@1.6) pkg> add GeophysicalModelGenerator
 ```
 You can test whether it works on your system with
-```julia
+```julia-repl
 julia> ]
-(@v1.10) pkg> test GeophysicalModelGenerator
+(@1.6) pkg> test GeophysicalModelGenerator
 ```
 and use it with
-```julia
+```julia-repl
 julia> using GeophysicalModelGenerator
 ```
 
diff --git a/docs/src/assets/img/Basic_Tutorial_Paraview_1.png b/docs/src/assets/img/Basic_Tutorial_Paraview_1.png
new file mode 100644
index 00000000..1b4cbd2f
Binary files /dev/null and b/docs/src/assets/img/Basic_Tutorial_Paraview_1.png differ
diff --git a/docs/src/assets/img/Basic_Tutorial_Paraview_2.png b/docs/src/assets/img/Basic_Tutorial_Paraview_2.png
new file mode 100644
index 00000000..fc25cb2a
Binary files /dev/null and b/docs/src/assets/img/Basic_Tutorial_Paraview_2.png differ
diff --git a/docs/src/assets/img/Basic_Tutorial_Paraview_3.png b/docs/src/assets/img/Basic_Tutorial_Paraview_3.png
new file mode 100644
index 00000000..79f5a940
Binary files /dev/null and b/docs/src/assets/img/Basic_Tutorial_Paraview_3.png differ
diff --git a/docs/src/man/Tutorial_Basic.md b/docs/src/man/Tutorial_Basic.md
index c01490f4..958e55c0 100644
--- a/docs/src/man/Tutorial_Basic.md
+++ b/docs/src/man/Tutorial_Basic.md
@@ -21,7 +21,7 @@ Tomo_Alps_full = load_GMG("https://zenodo.org/records/10738510/files/Paffrath_20
 ```
 
 ````
-GeoData 
+GeoData
   size      : (162, 130, 42)
   lon       ϵ [ -13.3019 : 35.3019]
   lat       ϵ [ 30.7638 : 61.2362]
@@ -49,7 +49,7 @@ Topo_Alps = load_GMG("https://zenodo.org/records/10738510/files/AlpsTopo.jld2?do
 ```
 
 ````
-GeoData 
+GeoData
   size      : (961, 841, 1)
   lon       ϵ [ 4.0 : 20.0]
   lat       ϵ [ 36.0 : 50.0]
@@ -70,9 +70,17 @@ Saved file: Topo_Alps.vts
 
 ````
 
-If we open both datasets in Paraview, we see this (after giving some color to the topography):
+If we open both datasets in Paraview, and changing both files from outline/solid colors to the corresponding data field, we see:
+![Basic_Tutorial_1](../assets/img/Basic_Tutorial_Paraview_1.png)
+Now we can change the colormap on the right side, marked by a red square. For topography we use the `Oleron` colormap, which you can download [here](https://www.fabiocrameri.ch/colourmaps/).
+For the tomography we use the `Roma` scientific colormap. You will now see a blue'ish box of the tomography, this is not the best color to visualise the data. Let's invert the colormap by clicking on the item marked by the blue arrow.
+Now we see the tomography in a more intuitive way, but the topography is not visible anymore. We can change the opacity of the tomography by setting a value in the `Opacity` field marked by the red square.
+Note that you will need to adapt the range of the topography colormap as the change in color is not at 0.0. By clicking on the item marked by the black arrow, you can set your desired range.
+
+![Basic_Tutorial_1](../assets/img/Basic_Tutorial_Paraview_2.png)
+
+Now you should see something like this:
 ![Basic_Tutorial_1](../assets/img/Basic_Tutorial_1.png)
-Note that I use the `Oleron` scientific colormap for the tomography which you can download [here](https://www.fabiocrameri.ch/colourmaps/)
 
 ### 2. Extract subset of data
 As you can see the tomographic data covers a much larger area than the Alps itself, and in most of that area there is no data.
@@ -89,7 +97,7 @@ Saved file: Tomo_Alps.vts
 
 ````
 
-Which looks like:
+After loading the new data again in paraview, switching to the proper data field and adjusting the colormap, you should see something like this:
 ![Basic_Tutorial_2](../assets/img/Basic_Tutorial_2.png)
 
 ### 3. Create cross sections
@@ -102,7 +110,7 @@ data_200km = cross_section(Tomo_Alps, Depth_level=-200)
 ```
 
 ````
-GeoData 
+GeoData
   size      : (54, 60, 1)
   lon       ϵ [ 3.9057 : 19.9057]
   lat       ϵ [ 35.9606 : 49.8976]
@@ -118,7 +126,7 @@ data_200km_exact = cross_section(Tomo_Alps, Depth_level=-200, Interpolate=true)
 ```
 
 ````
-GeoData 
+GeoData
   size      : (100, 100, 1)
   lon       ϵ [ 3.9057 : 19.9057]
   lat       ϵ [ 35.9606 : 49.8976]
@@ -128,7 +136,7 @@ GeoData
 ````
 
 In general, you can get help info for all functions with `?`:
-```julia
+```julia-repl
 help?> cross_section
 search: cross_section cross_section_volume cross_section_points cross_section_surface flatten_cross_section
 
@@ -176,7 +184,7 @@ Cross_vert = cross_section(Tomo_Alps, Start=(5,47), End=(15,44))
 ```
 
 ````
-GeoData 
+GeoData
   size      : (100, 100, 1)
   lon       ϵ [ 5.0 : 15.0]
   lat       ϵ [ 47.0 : 44.0]
@@ -198,8 +206,12 @@ Saved file: data_200km.vts
 
 ````
 
+After loading the data in Paraview, you can use the `Clip` tool on the topography to only show the topography above sealevel and make it 60% transparent. Also adjust the colormap of the tomography to 5.0 and -5.0
+
+![Basic_Tutorial_3](../assets/img/Basic_Tutorial_Paraview_3.png)
+
+After doing all these steps, you should see something like this:
 ![Basic_Tutorial_3](../assets/img/Basic_Tutorial_3.png)
-In creating this image, I used the `Clip` tool of Paraview to only show topography above sealevel and made it 50% transparent.
 
 ### 4. Cartesian data
 As you can see, the curvature or the Earth is taken into account here. Yet, for many applications it is more convenient to work in Cartesian coordinates (kilometers) rather then in geographic coordinates.
@@ -213,7 +225,7 @@ Topo_cart = convert2CartData(Topo_Alps, proj)
 ```
 
 ````
-CartData 
+CartData
     size    : (961, 841, 1)
     x       ϵ [ -748.7493528015041 : 695.3491277129204]
     y       ϵ [ -781.2344794653393 : 831.6826244089501]
@@ -229,7 +241,7 @@ Tomo_cart = convert2CartData(Tomo_Alps, proj)
 ```
 
 ````
-CartData 
+CartData
     size    : (54, 60, 39)
     x       ϵ [ -757.8031278236692 : 687.0608438357591]
     y       ϵ [ -785.601866956207 : 821.3433749818317]
@@ -269,7 +281,7 @@ Tomo_rect = project_CartData(Tomo_rect, Tomo_Alps, proj)
 ```
 
 ````
-CartData 
+CartData
     size    : (116, 121, 117)
     x       ϵ [ -550.0 : 600.0]
     y       ϵ [ -500.0 : 700.0]
@@ -286,7 +298,7 @@ Topo_rect = project_CartData(Topo_rect, Topo_Alps, proj)
 ```
 
 ````
-CartData 
+CartData
     size    : (1151, 1201, 1)
     x       ϵ [ -550.0 : 600.0]
     y       ϵ [ -500.0 : 700.0]
@@ -315,4 +327,3 @@ At this stage, the data can directly be used to generate cartesian numerical mod
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
-
diff --git a/docs/src/man/Tutorial_Chmy_MPI.md b/docs/src/man/Tutorial_Chmy_MPI.md
index fc078bfd..c4a57b87 100644
--- a/docs/src/man/Tutorial_Chmy_MPI.md
+++ b/docs/src/man/Tutorial_Chmy_MPI.md
@@ -156,9 +156,8 @@ You can than run it with:
 mpiexecjl -n 4 --project=. julia Tutorial_Chmy_MPI.jl
 ```
 
-The full file can be downloaded [here](../../../tutorials/Tutorial_Chmy_MPI.jl)
+The full file can be downloaded [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorials/Tutorial_Chmy_MPI.jl)
 
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
-
diff --git a/docs/src/man/installation.md b/docs/src/man/installation.md
index 38181b7d..b49afb13 100644
--- a/docs/src/man/installation.md
+++ b/docs/src/man/installation.md
@@ -15,7 +15,7 @@ You start julia on the command line with:
 kausb$ julia
 ```
 This will start the command-line interface of julia:
-```julia
+```julia-repl
                _
    _       _ _(_)_     |  Documentation: https://docs.julialang.org
   (_)     | (_) (_)    |
@@ -25,22 +25,22 @@ This will start the command-line interface of julia:
  _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
 |__/                   |
 
-julia> 
+julia>
 ```
 
 From the julia prompt, you start the package manager by typing `]`:
-```julia
-(@v1.6) pkg> 
+```julia-repl
+(@1.6) pkg>
 ```
 And you return to the command line with a backspace.
 
 Also useful is that julia has a build-in terminal, which you can reach by typing `;` on the command line:
-```julia
+```julia-repl
 julia>;
-shell> 
+shell>
 ```
 In the shell, you can use the normal commands like listing the content of a directory, or the current path:
-```julia
+```julia-repl
 shell> ls
 LICENSE         Manifest.toml   Project.toml    README.md       docs            src             test            tutorial
 shell> pwd
@@ -48,9 +48,9 @@ shell> pwd
 ```
 As before, return to the main command line (called `REPL`) with a backspace.
 
-If you want to see help information for any julia function, type `?` followed by the command. 
+If you want to see help information for any julia function, type `?` followed by the command.
 An example for `tan` is:
-```julia
+```julia-repl
 help?> tan
 search: tan tanh tand atan atanh atand instances transpose transcode contains UnitRange ReentrantLock StepRange StepRangeLen trailing_ones trailing_zeros
 
@@ -73,14 +73,14 @@ search: tan tanh tand atan atanh atand instances transpose transcode contains Un
   2×2 Matrix{Float64}:
    -1.09252  -1.09252
    -1.09252  -1.09252
-``` 
+```
 
 If you are in a directory that has a julia file (which have the extension `*.jl`), you can open that file with Visual Studio Code:
-```julia
+```julia-repl
 shell> code runtests.jl
 ```
 Execute the file with:
-```julia
+```julia-repl
 julia> include("runtests")
 ```
 Note that you do not include the `*.jl` extension.
@@ -88,44 +88,40 @@ Note that you do not include the `*.jl` extension.
 
 ### 4. Install GeophysicalModelGenerator.jl
 In order to install GeophysicalModelGenerator.jl, start julia and go to the package manager:
-```julia
+```julia-repl
 julia> ]
-(@v1.6) pkg> add GeophysicalModelGenerator
+(@v1.11) pkg> add GeophysicalModelGenerator
 ```
 This will automatically install various other packages it relies on (using the correct version).
 
 If you want, you can test if it works on your machine by running the test suite in the package manager:
-```julia
+```julia-repl
 julia> ]
-(@v1.6) pkg> test GeophysicalModelGenerator
+(@1.6) pkg> test GeophysicalModelGenerator
 ```
 Note that we run these tests automatically on Windows, Linux and Mac every time we add a new feature to GeophysicalModelGenerator (using different julia versions). This Continuous Integration (CI) ensures that new features do not break others in the package. The results can be seen [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/actions).
 
 The installation of `GMG` only needs to be done once, and will precompile the package and all other dependencies.
 
 If you, at a later stage, want to upgrade to the latest version of `GMG`, you can type:
-```julia
+```julia-repl
 julia> ]
-(@v1.6) pkg> update GeophysicalModelGenerator
+(@1.6) pkg> update GeophysicalModelGenerator
 ```
 
 You can load GeophysicalModelGenerator, for example to create cross-sections, with:
-```julia
+```julia-repl
 julia> using GeophysicalModelGenerator
 ```
 
 ### 5. Other useful packages
-As you will work your way through the tutorials you will see that we often use external packages, for example to load ascii data files into julia. You will find detailed instructions in the respective tutorials. 
+As you will work your way through the tutorials you will see that we often use external packages, for example to load ascii data files into julia. You will find detailed instructions in the respective tutorials.
 
 If you already want to install some of those, here our favorites. Install them through the package manager:
 
-- [CSV](https://github.com/JuliaData/CSV.jl): Read comma-separated data files into julia.  
-- [Plots](https://github.com/JuliaPlots/Plots.jl): Create all kinds of plots in julia (quite an extensive package, but very useful to have). 
+- [CSV](https://github.com/JuliaData/CSV.jl): Read comma-separated data files into julia.
+- [Plots](https://github.com/JuliaPlots/Plots.jl): Create all kinds of plots in julia (quite an extensive package, but very useful to have).
 - [JLD2](https://github.com/JuliaIO/JLD2.jl): This allows saving julia objects (such as a tomographic model) to a binary file and load it again at a later stage.
 - [Geodesy](https://github.com/JuliaGeo/Geodesy.jl): Convert UTM coordinates to latitude/longitude/altitude.
 - [NetCDF](https://github.com/JuliaGeo/NetCDF.jl): Read NetCDF files.
 - [GMT](https://github.com/GenericMappingTools/GMT.jl): A julia interface to the Generic Mapping Tools (GMT), which is a highly popular package to create (geophysical) maps. Note that installing `GMT.jl` is more complicated than installing the other packages listed above, as you first need to have a working version of `GMT` on your machine (it is not yet installed automatically). Installation instructions for Windows/Linux are on their webpage. On a mac, we made the best experiences by downloading the binaries from their webpage and not using a package manager to install GMT.
-
-
-
-
diff --git a/docs/src/man/projection.md b/docs/src/man/projection.md
index e8b4e9c9..f8934ffd 100644
--- a/docs/src/man/projection.md
+++ b/docs/src/man/projection.md
@@ -2,7 +2,7 @@
 
 Typically, you load a dataset by reading it into julia and either generating a `GeoData` structure (in case you have `longitude/latitude/depth` info), or as `UTMData` (in case the data is in `UTM coordinates`, which requires you to specify the zone & hemisphere).
 
-If you write the data to `Paraview`, it is internally converted to a Paraview structure (which involves `x,y,z` Cartesian Earth-Centered-Earth-Fixed (ECEF) coordinates using the `wgs84` ellipsoid). 
+If you write the data to `Paraview`, it is internally converted to a Paraview structure (which involves `x,y,z` Cartesian Earth-Centered-Earth-Fixed (ECEF) coordinates using the `wgs84` ellipsoid).
 
 Yet, if you do geodynamic calculations the chances are that the geodynamic code does not operate in spherical coordinates, but rather use cartesian ones. In that case you should transfer your data to the `CartData` structure, which requires you to specify a `ProjectionPoint` that is a point on the map that will later have the coordinates `(0,0)` in the `CartData` structure.
 
@@ -10,10 +10,10 @@ Yet, if you do geodynamic calculations the chances are that the geodynamic code
 #### 1. Converting
 Converting from one coordinate system to the other is straightforward. Let's use Europe as an example:
 
-```julia
+```julia-repl
 julia> using GeophysicalModelGenerator, GMT
 julia> Topo = import_topo(lon = [-10, 45], lat=[25, 50], file="@earth_relief_20m")
-GeoData 
+GeoData
   size  : (165, 75, 1)
   lon   ϵ [ -10.0 : 44.666666666666664]
   lat   ϵ [ 25.0 : 49.666666666666664]
@@ -22,13 +22,13 @@ GeoData
 julia> write_paraview(Topo,"Topo")
 Saved file: Topo.vts
 ```
-The result is shown on the globe as: 
+The result is shown on the globe as:
 ![Topo_Europe_GeoData](../assets/img/Topo_Europe_GeoData.png)
 
 You can convert this to UTM zone as:
-```julia
+```julia-repl
 julia> convert(UTMData, Topo)
-UTMData 
+UTMData
   UTM zone : 29-38 North
     size   : (165, 75, 1)
     EW     ϵ [ 197181.31221507967 : 769155.4572884373]
@@ -40,15 +40,15 @@ As the area is large, it covers a range of `UTM` zones (and every point has a UT
 
 Yet, what we could do instead is show all data with respect to a single UTM zone. For this, we have to select a point around which we project (in this case more or less in the center):
 
-```julia
+```julia-repl
 julia> p=ProjectionPoint(Lon=17.3, Lat=37.5)
 ProjectionPoint(37.5, 17.3, 703311.4380385976, 4.152826288024972e6, 33, true)
 ```
 
 Projecting the `GeoData` set using this projection point is done with:
-```julia
+```julia-repl
 julia> convert2UTMzone(Topo,p)
-UTMData 
+UTMData
   UTM zone : 33-33 North
     size   : (165, 75, 1)
     EW     ϵ [ -2.0750691599137965e6 : 3.581351293385453e6]
@@ -56,13 +56,13 @@ UTMData
     depth  ϵ [ -4985.5 m : 3123.0 m]
     fields : (:Topography,)
 ```
-Whereas this is now in UTM Data (in meters), it is distorted. 
+Whereas this is now in UTM Data (in meters), it is distorted.
 ![Topo_Europe_UTMData](../assets/img/Topo_Europe_UTMData.png)
 
 Often it is more convenient to have this in `CartData`, which is done in a similar manner:
-```julia
+```julia-repl
 julia> Topo_Cart = convert2CartData(Topo,p)
-CartData 
+CartData
     size   : (165, 75, 1)
     x      ϵ [ -2778.3805979523936 km : 2878.039855346856 km]
     y      ϵ [ -1387.8785405466067 km : 1785.2879241356998 km]
@@ -72,13 +72,13 @@ CartData
 This shows that the model is ~5600 by 3000 km.
 ![Topo_Europe_CartData](../assets/img/Topo_Europe_CartData.png)
 
-Whereas this is ok to look at and compare with a LaMEM model setup, we cannot use it to perform internal calculations (or to generate a LaMEM model setup), because the `x` and `y` coordinates are distorted and not orthogonal. 
+Whereas this is ok to look at and compare with a LaMEM model setup, we cannot use it to perform internal calculations (or to generate a LaMEM model setup), because the `x` and `y` coordinates are distorted and not orthogonal.
 
 #### 2. Projecting data
 For use with LaMEM, you would need an orthogonal cartesian grid. From the last command above we get some idea on the area, so we can create this:
-```julia
+```julia-repl
 julia> Topo_Cart_orth  = CartData(xyz_grid(-2000:20:2000,-1000:20:1000,0))
-CartData 
+CartData
     size   : (201, 101, 1)
     x      ϵ [ -2000.0 km : 2000.0 km]
     y      ϵ [ -1000.0 km : 1000.0 km]
@@ -86,20 +86,20 @@ CartData
     fields : (:Z,)
 ```
 Next, we can project the topographic data (in `GeoData` format) on this orthogonal grid
-```julia
+```julia-repl
 julia> Topo_Cart_orth  = project_CartData(Topo_Cart_orth, Topo, p)
-CartData 
+CartData
     size   : (201, 101, 1)
     x      ϵ [ -2000.0 km : 2000.0 km]
     y      ϵ [ -1000.0 km : 1000.0 km]
     z      ϵ [ -4.485650671162607 km : 2.5909655318121865 km]
     fields : (:Topography,)
-julia> write_paraview(Topo_Cart_orth,"Topo_Cart_orth");    
+julia> write_paraview(Topo_Cart_orth,"Topo_Cart_orth");
 ```
 ![Topo_Europe_CartData_Proj](../assets/img/Topo_Europe_CartData_Proj.png)
 So this interpolates the topographic data from the `GeoData` to the orthogonal cartesian grid (which can be used with LaMEM, for example).
 
-You can do similar projections with full 3D data sets or pointwise data. 
+You can do similar projections with full 3D data sets or pointwise data.
 
 #### 3. List of relevant functions
 
diff --git a/docs/src/man/tutorial_ISC_data.md b/docs/src/man/tutorial_ISC_data.md
index 998f34cc..17d39d66 100644
--- a/docs/src/man/tutorial_ISC_data.md
+++ b/docs/src/man/tutorial_ISC_data.md
@@ -10,7 +10,7 @@ You can get data from the ISC catalogue here:
 The catalogue will give you an on screen CSV output that will then have to be copied to a file of your choice (here we will call it `ISC1.dat`). Do that and start julia from the directory where it was downloaded.
 
 #### 2. Read data into Julia
-The main data-file, `ISC1.dat`, has 23 lines of comments (indicated with `#`), after which the data starts. We can use the julia package [https://github.com/JuliaData/CSV.jl](CSV.jl) to read in the data, and tell it that the data is separated by `,`.
+The main data-file, `ISC1.dat`, has 23 lines of comments (indicated with `#`), after which the data starts. We can use the julia package [CSV.jl](https://github.com/JuliaData/CSV.jl) to read in the data, and tell it that the data is separated by `,`.
 ```julia-repl
 julia> using CSV, GeophysicalModelGenerator
 julia> data_file = CSV.File("ISC1.dat",datarow=24,header=false,delim=',')
diff --git a/docs/src/man/tutorial_load3DSeismicData.md b/docs/src/man/tutorial_load3DSeismicData.md
index 7e054da1..12a5b0d4 100644
--- a/docs/src/man/tutorial_load3DSeismicData.md
+++ b/docs/src/man/tutorial_load3DSeismicData.md
@@ -2,18 +2,18 @@
 
 ## Goal
 This explains how to load a 3D P-wave model and plot it in Paraview as a 3D volumetric data set. It also shows how you can create horizontal or vertical cross-sections through the data in a straightforward manner and how you can extract subsets of the data;
-The example is the P-wave velocity model of the Alps as described in: 
+The example is the P-wave velocity model of the Alps as described in:
 
 Zhao, L., Paul, A., Malusà, M.G., Xu, X., Zheng, T., Solarino, S., Guillot, S., Schwartz, S., Dumont, T., Salimbeni, S., Aubert, C., Pondrelli, S., Wang, Q., Zhu, R., 2016. *Continuity of the Alpine slab unraveled by high-resolution P wave tomography*. Journal of Geophysical Research: Solid Earth 121, 8720–8737. [doi:10.1002/2016JB013310](https://doi.org/10.1002/2016JB013310)
 
 The data is given in ASCII format with longitude/latitude/depth/velocity anomaly (percentage) format.
 
 ## Steps
-#### 1. Download data 
+#### 1. Download data
 The data is can be downloaded from [https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/](https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/), where you should download the file `Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt`. Do that and start julia from the directory where it was downloaded.
 #### 2. Read data into Julia
 The dataset has no comments, and the data values in every row are separated by a space. In order to read this into julia as a matrix, we can use the build-in julia package `DelimitedFiles`.    We want the resulting data to be stored as double precision values (`Float64`), and the end of every line is a linebreak (`\n`).
-```julia
+```julia-repl
 julia> using DelimitedFiles
 julia> data=readdlm("Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt",' ',Float64,'\n', skipstart=0,header=false)
 1148774×4 Matrix{Float64}:
@@ -23,7 +23,7 @@ julia> data=readdlm("Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt",' ',Float64,'\n',
   0.45  38.0   -1001.0  -0.059
   0.6   38.0   -1001.0  -0.055
   0.75  38.0   -1001.0  -0.057
-  ⋮                     
+  ⋮
  17.25  51.95     -1.0  -0.01
  17.4   51.95     -1.0  -0.005
  17.55  51.95     -1.0   0.003
@@ -32,7 +32,7 @@ julia> data=readdlm("Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt",' ',Float64,'\n',
  18.0   51.95     -1.0   0.003
 ```
 Next, extract vectors from it:
-```julia
+```julia-repl
 julia> lon        = data[:,1];
 julia> lat        = data[:,2];
 julia> depth      = data[:,3];
@@ -43,7 +43,7 @@ Note that depth needs to with negative numbers.
 #### 3. Reformat the data
 
 Let's first have a look at the depth range of the data set:
-```julia
+```julia-repl
 julia> Depth_vec = unique(depth)
 101-element Vector{Float64}:
  -1001.0
@@ -61,20 +61,20 @@ julia> Depth_vec = unique(depth)
     -1.0
 ```
 So the data has a vertical spacing of 10 km.
-Next, let's check if the data is spaced in a regular manner in Lon/Lat direction. 
+Next, let's check if the data is spaced in a regular manner in Lon/Lat direction.
 For that, we read the data at a given depth level (say -101km) and plot it using the `Plots` package (you may have to install that first on your machine).
-```julia
+```julia-repl
 julia> using Plots
 julia> ind=findall(x -> x==-101.0, depth)
 julia> scatter(lon[ind],lat[ind],marker_z=dVp_perc[ind], ylabel="latitude",xlabel="longitude",markersize=2.5, c = :roma)
 ```
 ![DataPoints](../assets/img/Tutorial_Zhao_LatLon_data.png)
 
-Note that we employ the scientific colormap `roma` here. [This](https://docs.juliaplots.org/latest/generated/colorschemes/) gives an overview of available colormaps. You can download the colormaps for Paraview [here](https://www.fabiocrameri.ch/visualisation/).  
+Note that we employ the scientific colormap `roma` here. [This](https://docs.juliaplots.org/latest/generated/colorschemes/) gives an overview of available colormaps. You can download the colormaps for Paraview [here](https://www.fabiocrameri.ch/visualisation/).
 
 Clearly, the data is given as regular Lat/Lon points:
 
-```julia
+```julia-repl
 julia> unique(lon[ind])
 121-element Vector{Float64}:
   0.0
@@ -109,7 +109,7 @@ julia> unique(lat[ind])
 
 #### 3.1 Reshape data and save to paraview
 Next, we reshape the vectors with lon/lat/depth data into 3D matrixes:
-```julia
+```julia-repl
 julia> resolution =  (length(unique(lon)), length(unique(lat)), length(unique(depth)))
 (121, 94, 101)
 julia> Lon          = reshape(lon,      resolution);
@@ -119,7 +119,7 @@ julia> dVp_perc_3D  = reshape(dVp_perc, resolution);
 ```
 
 Check that the results are consistent
-```julia
+```julia-repl
 julia> iz=findall(x -> x==-101.0, Depth[1,1,:])
 1-element Vector{Int64}:
  91
@@ -130,10 +130,10 @@ heatmap(unique(lon), unique(lat),data[:,:,1]', c=:roma,title="$(Depth[1,1,iz]) k
 So this looks good.
 
 Next, we create a paraview file:
-```julia
+```julia-repl
 julia> using GeophysicalModelGenerator
 julia> Data_set    =   GeoData(Lon,Lat,Depth,(dVp_Percentage=dVp_perc_3D,))
-GeoData 
+GeoData
   size  : (121, 94, 101)
   lon   ϵ [ 0.0 - 18.0]
   lat   ϵ [ 38.0 - 51.95]
@@ -146,7 +146,7 @@ julia> write_paraview(Data_set, "Zhao_etal_2016_dVp_percentage")
 
 #### 4. Plotting data in Paraview
 In paraview, you should open the `*.vts` file, and press `Apply` (left menu) after doing that. Once you did that you can select `dVp_Percentage` and `Surface` (see red ellipses below).
-In paraview you can open the file and visualize it. 
+In paraview you can open the file and visualize it.
 ![Paraview_1](../assets/img/Tutorial_Zhao_Paraview_1.png)
 This visualisation employs the default colormap, which is not particularly good.
 
@@ -164,7 +164,7 @@ After pushing `Apply`, you'll see this:
 
 ![Paraview_5](../assets/img/Tutorial_Zhao_Paraview_5.png)
 
-If you want to plot iso-surfaces (e.g. at -3%), you can use the `Clip` option again, but this time select `scalar` and don't forget to unclick `invert`.
+If you want to plot iso-surfaces (e.g. at 3%), you can use the `Clip` option again, but this time select `scalar` and don't forget to unclick `invert`.
 
 ![Paraview_6](../assets/img/Tutorial_Zhao_Paraview_6.png)
 
@@ -175,9 +175,9 @@ In many cases you would like to create cross-sections through the 3D data sets a
 
 There is a simple way to achieve this using the `cross_section` function.
 To make a cross-section at a given depth:
-```julia
-julia> Data_cross  =   cross_section(Data_set, Depth_level=-100km)  
-GeoData 
+```julia-repl
+julia> Data_cross  =   cross_section(Data_set, Depth_level=-100km)
+GeoData
   size  : (121, 94, 1)
   lon   ϵ [ 0.0 : 18.0]
   lat   ϵ [ 38.0 : 51.95]
@@ -189,14 +189,14 @@ julia> write_paraview(Data_cross, "Zhao_CrossSection_100km")
 ```
 
 Or at a specific longitude:
-```julia
+```julia-repl
 julia> Data_cross  =   cross_section(Data_set, Lon_level=10)
-GeoData 
+GeoData
   size  : (1, 94, 101)
   lon   ϵ [ 10.05 : 10.05]
   lat   ϵ [ 38.0 : 51.95]
   depth ϵ [ -1001.0 km : -1.0 km]
-  fields: (:dVp_Percentage,) 
+  fields: (:dVp_Percentage,)
 julia> write_paraview(Data_cross, "Zhao_CrossSection_Lon10")
 1-element Vector{String}:
  "Zhao_CrossSection_Lon10.vts"
@@ -204,9 +204,9 @@ julia> write_paraview(Data_cross, "Zhao_CrossSection_Lon10")
 As you see, this cross-section is not taken at exactly 10 degrees longitude. That is because by default we don't interpolate the data, but rather use the closest point in longitude in the original data set.
 
 If you wish to interpolate the data, specify `Interpolate=true`:
-```julia
+```julia-repl
 julia> Data_cross = cross_section(Data_set, Lon_level=10, Interpolate=true)
-GeoData 
+GeoData
   size  : (1, 100, 100)
   lon   ϵ [ 10.0 : 10.0]
   lat   ϵ [ 38.0 : 51.95]
@@ -217,9 +217,9 @@ julia> write_paraview(Data_cross, "Zhao_CrossSection_Lon10_interpolated");
 as you see, this causes the data to be interpolated on a `(100,100)` grid (which can be changed by adding a `dims` input parameter).
 
 We can also create a diagonal cut through the model:
-```julia
+```julia-repl
 julia> Data_cross  =   cross_section(Data_set, Start=(1.0,39), End=(18,50))
-GeoData 
+GeoData
   size  : (100, 100, 1)
   lon   ϵ [ 1.0 : 18.0]
   lat   ϵ [ 39.0 : 50.0]
@@ -228,15 +228,15 @@ GeoData
 julia> write_paraview(Data_cross, "Zhao_CrossSection_diagonal")
 ```
 
-Here an image that shows the resulting cross-sections: 
+Here an image that shows the resulting cross-sections:
 ![Paraview_7](../assets/img/Tutorial_Zhao_Paraview_7.png)
 
 #### 6. Extract a (3D) subset of the data
 Sometimes, the data set covers a large region (e.g., the whole Earth), and you are only interested in a subset of this data for your project. You can obviously cut your data to the correct size in Paraview. Yet, an even easier way is the routine `extract_subvolume`:
 
-```julia
+```julia-repl
 julia> Data_subset     =   extract_subvolume(Data_set,Lon_level=(5,12), Lat_level=(40,45))
-GeoData 
+GeoData
   size  : (48, 35, 101)
   lon   ϵ [ 4.95 : 12.0]
   lat   ϵ [ 39.95 : 45.05]
@@ -249,9 +249,9 @@ This gives the resulting image. You can obviously use that new, smaller, data se
 
 By default, we extract the original data and do not interpolate it on a new grid.
 In some cases, you will want to interpolate the data on a different grid. Use the `Interpolate=true` option for that:
-```julia
+```julia-repl
 julia> Data_subset_interp     =   extract_subvolume(Data_set,Lon_level=(5,12), Lat_level=(40,45), Interpolate=true)
-GeoData 
+GeoData
   size  : (50, 50, 50)
   lon   ϵ [ 5.0 : 12.0]
   lat   ϵ [ 40.0 : 45.0]
@@ -261,21 +261,18 @@ julia> write_paraview(Data_subset, "Zhao_Subset_interp")
 ```
 #### 7. Load and save data to disk
 It would be useful to save the 3D data set we just created to disk, such that we can easily load it again at a later stage and create cross-sections etc, or compare it with other models. This can be done with the `save_GMG` function:
-```julia
+```julia-repl
 julia> save_GMG("Zhao_Pwave", Data_set)
 ```
 
 If you, at a later stage, want to load this file again do it as follows:
-```julia
+```julia-repl
 julia> Data_set_Zhao2016_Vp = load_GMG("Zhao_Pwave")
 ```
 
 
 #### 8. Julia script
 The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorial/Alps_VpModel_Zhao_etal_JGR2016.jl). You need to be in the same directory as in the data file, after which you can run it in julia with
-```julia
+```julia-repl
 include("Alps_VpModel_Zhao_etal_JGR2016.jl")
 ```
-
-
-
diff --git a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
index 3d3a4fa1..af8690bc 100644
--- a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
+++ b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
@@ -11,7 +11,7 @@ El-Sharkawy et al. (2020), *The Slab Puzzle of the Alpine‐Mediterranean Region
 The data is can be downloaded from [https://ds.iris.edu/files/products/emc/emc-files/El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc](https://ds.iris.edu/files/products/emc/emc-files/El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc). Do that and start julia from the directory where it was downloaded.
 
 #### 2. Read data into Julia
-The main data-file, `El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc`, is given as netCDF file. To read in data of this type, it is necessary to load an appropriate package. Here, we will use the [https://github.com/JuliaGeo/NetCDF.jl](NetCDF.jl) package. Download and install the package with:
+The main data-file, `El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc`, is given as netCDF file. To read in data of this type, it is necessary to load an appropriate package. Here, we will use the [NetCDF.jl](https://github.com/JuliaGeo/NetCDF.jl) package. Download and install the package with:
  ```julia
 julia> using Pkg
 julia> Pkg.add("NetCDF")
diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl
index b43c5b64..3a3ec732 100644
--- a/src/LaMEM_io.jl
+++ b/src/LaMEM_io.jl
@@ -97,12 +97,12 @@ Extracts a certain `keyword` from a LaMEM input `file` and convert it to a certa
 Optionally, you can also pass command-line arguments which will override the value read from the input file.
 
 # Example 1:
-```julia
+```julia-repl
 julia> nmark_z = ParseValue_LaMEM_InputFile("SaltModels.dat","nmark_z",Int64)
 ```
 
 # Example 2:
-```julia
+```julia-repl
 julia> nmark_z = ParseValue_LaMEM_InputFile("SaltModels.dat","nmark_z",Int64, args="-nel_x 128 -coord_x -4,4")
 ```
 
@@ -173,7 +173,7 @@ Parses a LaMEM input file and stores grid information in the `Grid` structure.
 Optionally, you can pass LaMEM command-line arguments as well.
 
 # Example 1
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("SaltModels.dat")
 LaMEM Grid:
 nel         : (32, 32, 32)
@@ -185,7 +185,7 @@ z           ϵ [-2.0 : 0.0]
 ```
 
 # Example 2 (with command-line arguments)
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("SaltModels.dat", args="-nel_x 64 -coord_x -4,4")
 LaMEM Grid:
   nel         : (64, 32, 32)
@@ -823,7 +823,7 @@ end
 Reads a parallel, rectilinear, `*.vts` file with the name `fname` and located in `dir` and create a 3D `Data` struct from it.
 
 # Example
-```julia
+```julia-repl
 julia> Data = read_data_PVTR("Haaksbergen.pvtr", "./Timestep_00000005_3.35780500e-01/")
 ParaviewData
   size  : (33, 33, 33)
@@ -1025,6 +1025,6 @@ function coordinate_grids(Data::LaMEM_grid; cell=false)
     if cell
         X,Y,Z = average_q1(X),average_q1(Y), average_q1(Z)
     end
-    
+
     return X,Y,Z
 end
diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index ba29e034..1771c60c 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -16,7 +16,7 @@ export  add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_po
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp, LinearWeightedTemperature,
         McKenzie_subducting_slab,
         ConstantPhase, LithosphericPhases,
-        Trench,
+        Trench, compute_slab_surface,
         compute_thermal_structure, compute_phase
 
 """
@@ -71,7 +71,7 @@ Example
 ========
 
 Example: Box with striped phase and constant temperature & a dip angle of 10 degrees:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -175,7 +175,7 @@ Examples
 ========
 
 Example 1) Box with constant phase and temperature & a dip angle of 10 degrees:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -194,7 +194,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 ```
 
 Example 2) Box with halfspace cooling profile
-```julia
+```julia-repl
 julia> Grid = CartData(xyz_grid(-1000:10:1000,0,-660:10:0))
 julia> Phases = zeros(Int32,   size(Grid));
 julia> Temp   = zeros(Float64, size(Grid));
@@ -291,7 +291,7 @@ Examples
 ========
 
 Example 1) Layer with constant phase and temperature
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -310,7 +310,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to para
 ```
 
 Example 2) Box with halfspace cooling profile
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 julia> Phases = zeros(Int32,   size(Grid.X));
 julia> Temp   = zeros(Float64, size(Grid.X));
@@ -393,7 +393,7 @@ Example
 ========
 
 Sphere with constant phase and temperature:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -462,7 +462,7 @@ Example
 ========
 
 Ellipsoid with constant phase and temperature, rotated 90 degrees and tilted by 45 degrees:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -534,7 +534,7 @@ Parameters
 ====
 - `Phase` - Phase array (consistent with Grid)
 - `Temp`  - Temperature array (consistent with Grid)
-- `Grid` - Grid structure (usually obtained with read_LaMEM_inputfile)
+- `Grid` - Grid structure (usually obtained with `read_LaMEM_inputfile`)
 - `base` - center coordinate of bottom of cylinder
 - `cap` - center coordinate of top of cylinder
 - `radius` - radius of the cylinder
@@ -547,7 +547,7 @@ Example
 ========
 
 Cylinder with constant phase and temperature:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -644,7 +644,7 @@ Example
 
 Polygon with constant phase and temperature:
 
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -825,7 +825,7 @@ Example
 ========
 
 Cylinder with constant phase and temperature:
-```julia
+```julia-repl
 julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
 LaMEM Grid:
   nel         : (32, 32, 32)
@@ -1431,12 +1431,12 @@ Thermal structure by McKenzie for a subducted slab that is fully embedded in the
 
 Parameters
 ===
-- Tsurface:     Top T [C]
-- Tmantle:      Bottom T [C]
-- Adiabat:      Adiabatic gradient in K/km
-- v_cm_yr:      Subduction velocity [cm/yr]
-- κ:            Thermal diffusivity [m2/s]
-- it:           Number iterations employed in the harmonic summation
+- `Tsurface`:     Top T [C]
+- `Tmantle`:      Bottom T [C]
+- `Adiabat`:      Adiabatic gradient in K/km
+- `v_cm_yr`:      Subduction velocity [cm/yr]
+- `κ`:            Thermal diffusivity [m2/s]
+- `it`:           Number iterations employed in the harmonic summation
 
 """
 @with_kw_noshow mutable struct McKenzie_subducting_slab <: AbstractThermalStructure
@@ -1458,12 +1458,12 @@ that the bottom of the slab is the coordinate Z=0. Internally the function shift
 Parameters
 
 =============================
-Temp Temperature array
-- `X`    X Array
-- `Y`    Y Array
-- `Z`    Z Array
-- `Phase` Phase array
-- `s`    McKenzie_subducting_slab
+- `Temp`:  Temperature array
+- `X`:    X Array
+- `Y`:    Y Array
+- `Z`:    Z Array
+- `Phase`: Phase array
+- `s`:    `McKenzie_subducting_slab`
 """
 function compute_thermal_structure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab)
     @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
@@ -1525,15 +1525,19 @@ end
     compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)
 
 Weight average along distance
+
 Do a weight average between two field along a specified direction
-Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of
-F1, while F2 decreases.
-This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in
+
+Given a distance (could be any array, from X,Y) -> the weight of F1 increase from the origin, while F2 decreases.
+
+This function has been conceived for averaging the solution of McKenzie and half space cooling models, but it
 can be used to smooth the temperature field from continent ocean:
--> Select the boundary to apply;
--> transform the coordinate such that dist represent the perpendicular direction along which you want to apply
-this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0;
--> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F.
+- Select the boundary to apply;
+- transform the coordinate such that dist represent the perpendicular direction along which you want to apply this smoothening
+  and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0;
+- Select the points that belongs to this area
+- compute the thermal fields {F1} {F2}
+- then modify F.
 """
 function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)
     @unpack w_min, w_max, crit_dist,dir = s;
@@ -1864,7 +1868,7 @@ Examples
 ========
 
 Example 1) Slab
-```julia
+```julia-repl
 julia> x     = LinRange(0.0,1200.0,128);
 julia> y     = LinRange(0.0,1200.0,128);
 julia> z     = LinRange(-660,50,128);
diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl
index a2af9da9..20633cd0 100644
--- a/src/WaterFlow.jl
+++ b/src/WaterFlow.jl
@@ -19,7 +19,7 @@ function spacing(lon,lat)
      dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2
      dlat[:,1] = dlat[:,2]
      dlat[:,end] = dlat[:,end-1]
-    
+
     return dlon, dlat
 end
 
@@ -39,26 +39,26 @@ end
 
 
 """
-    Topo_water, sinks, pits, bnds  = waterflows(Topo::GeoData; 
+    Topo_water, sinks, pits, bnds  = waterflows(Topo::GeoData;
         flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true,
         rainfall = nothing,
         minsize=300)
 
 Takes a GMG GeoData object of a topographic map and routes water through the grid. Optionally,
-you can specify `rainfall` in which case we accumulate the rain as specified in this 2D array instead of the cellarea. 
+you can specify `rainfall` in which case we accumulate the rain as specified in this 2D array instead of the cellarea.
 This allows you to, for example, sum, up water if you have variable rainfall in the area.
 The other options are as in the `waterflows` function of the package `WhereTheWaterFlows`.
 
 Example
 ===
-```julia
+```julia-repl
 # Download some topographic data
 julia> Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
 
 # Flow the water through the area:
 julia> Topo_water, sinks, pits, bnds  = waterflows(Topo)
 julia> Topo_water
-GeoData 
+GeoData
   size      : (961, 481, 1)
   lon       ϵ [ 6.5 : 7.3]
   lat       ϵ [ 50.2 : 50.59999999999999]
@@ -68,7 +68,7 @@ GeoData
 ```
 
 """
-function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing, minsize=300) 
+function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing, minsize=300)
 
     cellarea = cell_area(Topo)
     cellarea_m2 = cellarea
@@ -96,10 +96,10 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature;
     largest_catchment = catchment_large .== catchment_large[id_max]
     largest_area = copy(area)
     largest_area[.!largest_catchment] .= NaN
-    
+
     log10_area = log10.(area)
     log10_largest_area = log10.(largest_area)
 
-    Topo_water =  addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area, log10_largest_area))                        
-    return Topo_water, sinks, pits, bnds 
-end
\ No newline at end of file
+    Topo_water =  addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area, log10_largest_area))
+    return Topo_water, sinks, pits, bnds
+end
diff --git a/src/surface_functions.jl b/src/surface_functions.jl
index 17265abc..f69fc452 100644
--- a/src/surface_functions.jl
+++ b/src/surface_functions.jl
@@ -17,12 +17,12 @@ function is_surface(surf::AbstractGeneralGrid)
     return issurf
 end
 
-function  +(a::_T, b::_T) where _T<:AbstractGeneralGrid 
+function  +(a::_T, b::_T) where _T<:AbstractGeneralGrid
     @assert size(a) == size(b)
     return _addSurfaces(a,b)
 end
 
-function  -(a::_T, b::_T) where _T<:AbstractGeneralGrid 
+function  -(a::_T, b::_T) where _T<:AbstractGeneralGrid
     @assert size(a) == size(b)
     return _subtractSurfaces(a,b)
 end
@@ -71,11 +71,11 @@ This drapes fields of a data set `Data` on the topography `Topo`
 function drape_on_topo(Topo::GeoData, Data::GeoData)
     @assert is_surface(Topo)
     @assert is_surface(Data)
-   
+
     Lon,Lat,_   =   lonlatdepth_grid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]);
 
     # use nearest neighbour to interpolate data
-    idx         =   nearest_point_indices(Lon,Lat, vec(Data.lon.val), vec(Data.lat.val) ); 
+    idx         =   nearest_point_indices(Lon,Lat, vec(Data.lon.val), vec(Data.lat.val) );
 
     idx_out     =   findall(  (Lon .<  minimum(Data.lon.val)) .| (Lon .>  maximum(Data.lon.val)) .|
                               (Lat .<  minimum(Data.lat.val)) .| (Lat .>  maximum(Data.lat.val)) )
@@ -137,7 +137,7 @@ Drapes Cartesian Data on topography
 function drape_on_topo(Topo::CartData, Data::CartData)
     @assert is_surface(Topo)
     @assert is_surface(Data)
-    
+
     Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields )
     Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields )
 
@@ -152,12 +152,12 @@ end
 """
     surf_new = fit_surface_to_points(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)
 
-This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`) 
+This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`)
 
 """
 function fit_surface_to_points(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)
     @assert is_surface(surf)
-    
+
     idx = nearest_point_indices(NumValue(surf.lon),NumValue(surf.lat),  lon_pt, lat_pt);
     depth = NumValue(surf.depth)
     depth[idx] .= depth_pt[idx];
@@ -171,12 +171,12 @@ end
 """
     surf_new = fit_surface_to_points(surf::CartData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)
 
-This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`) 
+This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`)
 
 """
 function fit_surface_to_points(surf::CartData, X_pt::Vector, Y_pt::Vector, Z_pt::Vector)
     @assert is_surface(surf)
-    
+
     idx = nearest_point_indices(NumValue(surf.x),NumValue(surf.y),  X_pt[:], Y_pt[:]);
     depth = NumValue(surf.z)
     depth = Z_pt[idx]
@@ -197,7 +197,7 @@ This can be used, for example, to mask points above/below the Moho in a volumetr
 
 # Example
 First we create a 3D data set and a 2D surface:
-```julia
+```julia-repl
 julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
 julia> Data            =   Depth*2;
 julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon))
@@ -217,7 +217,7 @@ julia> Data_Moho       =   GeoData(Lon,Lat,Depth+Lon*km, (MohoDepth=Depth,))
     fields: (:MohoDepth,)
 ```
 Next, we intersect the surface with the data set:
-```julia
+```julia-repl
 julia> Above       =   above_surface(Data_set3D, Data_Moho);
 ```
 Now, `Above` is a boolean array that is true for points above the surface and false for points below and at the surface.
@@ -337,7 +337,7 @@ end
 
 Interpolates a 3D data set `V` on a surface defined by `Surf`.
 # Example
-```julia
+```julia-repl
 julia> Data
 ParaviewData
   size  : (33, 33, 33)
@@ -403,4 +403,4 @@ function interpolate_data_surface(V::GeoData, Surf::GeoData)
     Surf_interp = interpolate_datafields(V, Surf.lon.val, Surf.lat.val, Surf.depth.val)
 
     return Surf_interp
-end
\ No newline at end of file
+end
diff --git a/tutorials/Tutorial_AlpineData.jl b/tutorials/Tutorial_AlpineData.jl
index ea450cba..3f21a4da 100644
--- a/tutorials/Tutorial_AlpineData.jl
+++ b/tutorials/Tutorial_AlpineData.jl
@@ -10,9 +10,9 @@
 # This is a rather lengthy tutorial that combines different other tutorials, but it will guide you through all the steps necessary to obtain a somewhat comprehensive view of the European Alps and their subsurface from a geodynamical point of view.
 
 # ## 1. Surface Topography
-# In many cases, we want to add topographic data to our visualization. Here we use [GMT.jl](https://github.com/GenericMappingTools/GMT.jl) to download data from a certain region, and transfer that to `GMG`. 
+# In many cases, we want to add topographic data to our visualization. Here we use [GMT.jl](https://github.com/GenericMappingTools/GMT.jl) to download data from a certain region, and transfer that to `GMG`.
 # To add the GMT package, simply add it with the julia package manager:
-# ```julia
+# ```julia-repl
 # julia> ]
 # (@v1.10) pkg> add GMT
 # ```
@@ -23,25 +23,25 @@ using GeophysicalModelGenerator, GMT
 # When loading both packages, several `GMT` routines within `GMG` will be loaded. One of these routines is the function `import_topo`, where one simply has to provide the region for which to download the topographic data and the data source.
 
 Topo = import_topo([4,20,37,50], file="@earth_relief_01m.grd")
-# The data is available in different resolutions; see [here](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html) for an overview. Generally, it is advisable to not use the largest resolution if you have a large area, as the files become very large. 
+# The data is available in different resolutions; see [here](http://gmt.soest.hawaii.edu/doc/latest/grdimage.html) for an overview. Generally, it is advisable to not use the largest resolution if you have a large area, as the files become very large.
 
-# If you have issues with loading the topography with `GMT`, there is also the alternative to download the data yourself and import it using `Rasters.jl`. 
+# If you have issues with loading the topography with `GMT`, there is also the alternative to download the data yourself and import it using `Rasters.jl`.
 
 # We can now export this data to a `VTK` format so that we can visualize it with `Paraview`. To do so, `GMG` provides the function `write_paraview`:
-write_paraview(Topo, "Topography_Alps") 
+write_paraview(Topo, "Topography_Alps")
 
 # Also, if you want to save this data for later use in julia, you can save it as `*.jld2` file using the function `save_GMG`:
 save_GMG("Topography_Alps",Topo)
 
 # The result looks like:
-# ![Alps_Tutorial_1](../assets/img/Tut_Alp_Image1.png) 
+# ![Alps_Tutorial_1](../assets/img/Tut_Alp_Image1.png)
 # Note that we used the `Oleron` scientific colormap from [here](https://www.fabiocrameri.ch/colourmaps.php) for the topography.
 
 
 # ##  2. Moho topography
-# When looking at data concerning the Alpine subsurface, we are often interested in the depth of the [Moho](https://en.wikipedia.org/wiki/Mohorovi%C4%8Di%C4%87_discontinuity). 
+# When looking at data concerning the Alpine subsurface, we are often interested in the depth of the [Moho](https://en.wikipedia.org/wiki/Mohorovi%C4%8Di%C4%87_discontinuity).
 # ### 2.1 Download and import of the data
-# Here, we will use the dataset from Mroczek et al. (2023). This dataset is publicly available and can be downloaded from [here](https://datapub.gfz-potsdam.de/download/10.5880.GFZ.2.4.2021.009NUEfb/2021-009_Mroczek-et-al_SWATHD_moho_jul22.csv). 
+# Here, we will use the dataset from Mroczek et al. (2023). This dataset is publicly available and can be downloaded from [here](https://datapub.gfz-potsdam.de/download/10.5880.GFZ.2.4.2021.009NUEfb/2021-009_Mroczek-et-al_SWATHD_moho_jul22.csv).
 # To allow for downloading such data, we use the julia package `Downloads.jl`, which is a dependency of `GMG`. To download the Moho data in the current directory, simply type:
 download_data("https://datapub.gfz-potsdam.de/download/10.5880.GFZ.2.4.2021.009NUEfb/2021-009_Mroczek-et-al_SWATHD_moho_jul22.csv","MohoMroczek2023.csv")
 
@@ -51,7 +51,7 @@ using DelimitedFiles
 data_mroczek = readdlm("MohoMroczek2023.csv",',',header=false,skipstart=11)
 # Note that we skipped the first 11 lines of the file as they contain the file header. The result of this operation will look like this:
 #
-# ```julia
+# ```julia-repl
 # julia> data_mroczek = readdlm("MohoMroczek2023.csv",',',header=false,skipstart=11)
 # 40786×10 Matrix{Any}:
 #      "X"      "Y"      "Z"    "lat"    "lon"    "depth"   "tPs"   "k"   "interp"  "tag"
@@ -69,7 +69,7 @@ data_mroczek = readdlm("MohoMroczek2023.csv",',',header=false,skipstart=11)
 #  4186.85  1008.02  4639.51  47.1319  13.5368  40.8435    4.75    1.63  0          "PA"
 #  4188.49  1008.31  4638.04  47.1118  13.5355  40.7913    4.74    1.63  0          "PA"
 #  4190.14  1008.6   4636.57  47.0917  13.5341  40.7305    4.73    1.63  0          "PA"
-#     ⋮                                          ⋮                                  
+#     ⋮                                          ⋮
 #  4123.28  1111.52  4669.18  47.5537  15.0867  43.4301    5       1.67  0          "EU"
 #  4124.02  1111.57  4666.69  47.5336  15.0848  44.7763    5.13    1.67  0          "EU"
 #  4124.67  1111.59  4664.11  47.5136  15.0828  46.2535    5.28    1.67  0          "EU"
@@ -86,7 +86,7 @@ data_mroczek = readdlm("MohoMroczek2023.csv",',',header=false,skipstart=11)
 #  4126.5   1113.11  4643.52  47.3729  15.096   59.9569    6.71    1.66  0          "EU"
 # ```
 
-# We are now only interested in the depth of the Moho at a given longitude/latitude. To obtain these values, we now have to extract columns 4-6. In addition, we also extract the 10th column, as it contains an identifier for the tectonic unit the respective point belongs to. 
+# We are now only interested in the depth of the Moho at a given longitude/latitude. To obtain these values, we now have to extract columns 4-6. In addition, we also extract the 10th column, as it contains an identifier for the tectonic unit the respective point belongs to.
 lon        = zeros(size(data_mroczek,1)-1);     lon   .= data_mroczek[2:end,5];
 lat        = zeros(size(data_mroczek,1)-1);     lat   .= data_mroczek[2:end,4];
 depth      = zeros(size(data_mroczek,1)-1);     depth .= -1.0*data_mroczek[2:end,6]; #multiplied with -1, as we consider depth to be negative
@@ -107,18 +107,18 @@ using NearestNeighbors
 # Now that we have generated the grid, we can loop over our different tectonic units, extract the relevant data points and interpolate them to the regular grid:
 for iunit = 1:length(units)
     Dist              = zeros(size(Lon))
-    
+
     #Get all points belonging to the unit
     ind_unit    = findall( x -> occursin(units[iunit], x), tag) #index of the points belonging to that unit
     lon_tmp     = lon[ind_unit]
     lat_tmp     = lat[ind_unit]
     depth_tmp   = depth[ind_unit]
-    
-    #for later checking, we can now save the original point data as a VTK file: 
+
+    #for later checking, we can now save the original point data as a VTK file:
     data_Moho = GeophysicalModelGenerator.GeoData(lon_tmp,lat_tmp,depth_tmp,(MohoDepth=depth_tmp*km,))
     filename = "Mroczek_Moho_" * units[iunit]
     write_paraview(data_Moho, filename, PointsData=true)
-    
+
     #Now we create a KDTree for an effective nearest neighbor determination;
     kdtree = KDTree([lon_tmp';lat_tmp']; leafsize = 10)
     points = [Lon[:]';Lat[:]']
@@ -126,17 +126,17 @@ for iunit = 1:length(units)
     dists = reduce(vcat,dists)
     idxs  = reduce(vcat,idxs)
     idxs  = reduce(vcat,idxs)
-    
+
     #Having determined the nearest neighbor for each point in our regular grid, we can now directly assign the respective depth. Whenever the nearest neighbor is further than a certain distance away, we assume that there is no Moho at this point and do not assign a depth to that point.
     for i=1:length(idxs)
-        if dists[i]<0.02 
+        if dists[i]<0.02
             Depth[i] = depth_tmp[idxs[i]]*km
         else
             Depth[i] = NaN*km
         end
-        Dist[i]  = dists[i]    
+        Dist[i]  = dists[i]
     end
-    
+
     #As we will be using the data later, we would also like to provide some Metadata so that we know where it is coming from:
     Data_attribs   = Dict(
         "author"=>  "Mroczek et al.",
@@ -144,25 +144,25 @@ for iunit = 1:length(units)
         "doi"=>"https://doi.org/10.5880/GFZ.2.4.2021.009",
         "url"=>"https://nextcloud.gfz-potsdam.de/s/zB5dPNby6X2Kjnj",
     )
-    
+
     #Finally, we can now export that data to VTK and save a `jld2` file using the `save_GMG` routine
     Data_Moho = GeophysicalModelGenerator.GeoData(Lon, Lat, Depth, (MohoDepth=Depth,PointDist=Dist),Data_attribs)
     filename = "Mrozek_Moho_Grid_" * units[iunit]
     write_paraview(Data_Moho, filename)
     save_GMG(filename,Topo)
-    
+
 end
-    
+
 # Just for checking, we can also plot both the original data and the resulting interpolated Moho:
-# ![Alps_Tutorial_2](../assets/img/Tut_Alp_Image2.png) 
+# ![Alps_Tutorial_2](../assets/img/Tut_Alp_Image2.png)
 
 
 # ## 3. Seismicity
-# Earthquakes are always interesting, so lets import the seismicity data from ISC. 
+# Earthquakes are always interesting, so lets import the seismicity data from ISC.
 # ### 3.1 Download and import
-# ISC provides a method to download parts of it's catalogue via a web interface. 
+# ISC provides a method to download parts of it's catalogue via a web interface.
 # See the description of the interface [here](http://www.isc.ac.uk/iscbulletin/search/webservices/bulletin/).
-# We will now download all reviewed earthquake data between 1990 and 2015 in the same region as the extracted topography. 
+# We will now download all reviewed earthquake data between 1990 and 2015 in the same region as the extracted topography.
 # We will only consider earthquakes with a magnitude larger than 3. The resulting dataset is quite large, so consider to either limit the time range or the magnitude range.
 download_data("http://www.isc.ac.uk/cgi-bin/web-db-run?request=COLLECTED&req_agcy=ISC-EHB&out_format=QuakeML&ctr_lat=&ctr_lon=&radius=&max_dist_units=deg&searchshape=RECT&top_lat=49&bot_lat=37&left_lon=4&right_lon=20&srn=&grn=&start_year=1990&start_month=1&start_day=01&start_time=00%3A00%3A00&end_year=2015&end_month=12&end_day=31&end_time=00%3A00%3A00&min_dep=&max_dep=&min_mag=3.0&max_mag=&req_mag_type=Any&req_mag_agcy=Any&min_def=&max_def=&prime_only=on&include_magnitudes=on&table_owner=iscehb","ISCData.xml")
 
@@ -172,15 +172,15 @@ Data_ISC = getlonlatdepthmag_QuakeML("ISCData.xml");
 write_paraview(Data_ISC, "EQ_ISC", PointsData=true);
 save_GMG("EQ_ISC",Data_ISC)
 
-# ![Alps_Tutorial_3](../assets/img/Tut_Alp_Image3.png) 
+# ![Alps_Tutorial_3](../assets/img/Tut_Alp_Image3.png)
 
 
 # ## 4. GPS data
-# Besides data on the structure of the subsurface, it is also nice to see the dynamics of a region. Dynamic processes can be nicely seen in the surface velocities given by GPS data. 
+# Besides data on the structure of the subsurface, it is also nice to see the dynamics of a region. Dynamic processes can be nicely seen in the surface velocities given by GPS data.
 # As GPS data consists of three-dimensional vectors, we have to treat it differently than the seismicity data in the previous section. The example is based on a paper by Sanchez et al. (2018) [https://essd.copernicus.org/articles/10/1503/2018/#section7](https://essd.copernicus.org/articles/10/1503/2018/#section7).
 #
-# ### 4.1. Download and import GPS data: 
-# The data related to the paper can be downloaded from: [here](https://doi.pangaea.de/10.1594/PANGAEA.886889). There you will find links to several data sets. 
+# ### 4.1. Download and import GPS data:
+# The data related to the paper can be downloaded from: [here](https://doi.pangaea.de/10.1594/PANGAEA.886889). There you will find links to several data sets.
 # Some are the data on the actual stations and some are interpolated data on a grid. Here, we will use the gridded data as an example, and will therefore download the following data sets:
 #
 # - ALPS2017_DEF_HZ	Surface deformation model of the Alpine Region	[https://store.pangaea.de/Publications/Sanchez-etal_2018/ALPS2017_DEF_HZ.GRD](https://store.pangaea.de/Publications/Sanchez-etal_2018/ALPS2017_DEF_HZ.GRD)
@@ -214,7 +214,7 @@ Plots.scatter(lon_vz,lat_vz)
 nlon = length(unique(lon_vz))
 nlat = length(unique(lat_vz))
 
-# So we have a `41` by `31` grid. GMG requires 3D matrixes for the data (as we want to plot the results in Paraview in 3D). 
+# So we have a `41` by `31` grid. GMG requires 3D matrixes for the data (as we want to plot the results in Paraview in 3D).
 # That is why we first initialize 3D matrixes for `lon,lat,Vz`:
 Lon   =   zeros(nlon,nlat,1)
 Lat   =   zeros(nlon,nlat,1)
@@ -244,10 +244,10 @@ vn     =   data_vh[:,4]
 
 # Let's have a look at how the data points of this dataset are distributed:
 Plots.scatter(lon_vh,lat_vh)
-# So it appears that the horizontal velocities are given on the same regular grid as well, 
-# but not in regions which are covered with water. 
-# This thus requires a bit more work to transfer them to a rectangular grid. 
-# The strategy we take is to first define 2D matrixes with horizontal velocities with the same size as `Vz` 
+# So it appears that the horizontal velocities are given on the same regular grid as well,
+# but not in regions which are covered with water.
+# This thus requires a bit more work to transfer them to a rectangular grid.
+# The strategy we take is to first define 2D matrixes with horizontal velocities with the same size as `Vz`
 # initialized with `NaN` (not a number).
 Ve = fill(NaN,size(Vz));
 Vn = fill(NaN,size(Vz));
@@ -258,18 +258,18 @@ for i in eachindex(lon_vh)
     Vn[ind] .= vn[i];
 end
 
-# At this stage, we have horizontal and vertical velocities in units of `m/yr`. 
+# At this stage, we have horizontal and vertical velocities in units of `m/yr`.
 # Yet, given the small velocities in the Alps, it makes more sense to have them in units of `mm/yr`:
 Vz = Vz*1000;
 Ve = Ve*1000;
 Vn = Vn*1000;
 
 # And their magnitude is
-Vmagnitude  =   sqrt.(Ve.^2 + Vn.^2 + Vz.^2); 
+Vmagnitude  =   sqrt.(Ve.^2 + Vn.^2 + Vz.^2);
 
 # ### 4.2 Interpolate topography on the grid
-# At this stage we have the 3D velocity components on a grid. Yet, we don't have information yet about the elevation of the stations (as the provided data set did not give this). 
-# We could ignore that and set the elevation to zero, which would allow saving the data directly. 
+# At this stage we have the 3D velocity components on a grid. Yet, we don't have information yet about the elevation of the stations (as the provided data set did not give this).
+# We could ignore that and set the elevation to zero, which would allow saving the data directly.
 # Yet, a better way is to load the topographic map of the area and interpolate the elevation to the velocity grid.
 # As we have already the loaded the topographic map in section 1 of this tutorial, we can simply reuse it. To interpolate, we will use the function `interpolate_datafields_2D`
 
@@ -285,11 +285,11 @@ Data_GPS_Sanchez = GeoData(Lon,Lat,topo_v,(Velocity_mm_year=(Ve,Vn,Vz),V_north=V
 write_paraview(Data_GPS_Sanchez, "GPS_Sanchez")
 save_GMG("GPS_Sanchez",Data_GPS_Sanchez)
 
-# ![Alps_Tutorial_4](../assets/img/Tut_Alp_Image4.png) 
+# ![Alps_Tutorial_4](../assets/img/Tut_Alp_Image4.png)
 
 # ## 5. Seismic tomography data
 # Finally, we'd like to have a look at the subsurface by looking at a seismic tomography. To do so, we'll first download the tomography published by [Rappisi et al.(2022)](https://doi.org/10.1029/2021JB023488).
-# The data is provided as `NetCDF` files: 
+# The data is provided as `NetCDF` files:
 download_data("https://figshare.com/ndownloader/files/34093955","aniNEWTON21.nc")
 
 # We can load this file with the `NCDatasets` package.
@@ -303,7 +303,7 @@ dlnVp   = dataset["dlnVp"]
 Vp      = dataset["Vpi"]
 Depth   = dataset["Zg"]
 
-# As longitude and latitude are only given as 2D grids, we here have to convert them to 3D matrices. 
+# As longitude and latitude are only given as 2D grids, we here have to convert them to 3D matrices.
 Lon = repeat(lon[:,:],1,1,size(Depth,3));
 Lat = repeat(lat[:,:],1,1,size(Depth,3));
 
@@ -326,13 +326,13 @@ write_paraview(Data, "Rappisi2022")
 save_GMG("Rappisi2022",Data)
 
 # The result looks like:
-# ![Alps_Tutorial_5](../assets/img/Tut_Alp_Image5.png) 
+# ![Alps_Tutorial_5](../assets/img/Tut_Alp_Image5.png)
 
 # For the sake of this tutorial, we have now imported all the data we would like to look at. All that is missing is now a joint visualization
-# of these datasets. To obtain this visualization, we will load all the `VTK` files into Paraview and have a look: 
-# ![Alps_Tutorial_6](../assets/img/GMG_AlpineData.png) 
+# of these datasets. To obtain this visualization, we will load all the `VTK` files into Paraview and have a look:
+# ![Alps_Tutorial_6](../assets/img/GMG_AlpineData.png)
 #
-# A Paraview statefile that reproduces this visualization is available under `tutorials/Tutorial_AlpineData.pvsm`. 
+# A Paraview statefile that reproduces this visualization is available under `tutorials/Tutorial_AlpineData.pvsm`.
 
 #src Note: The markdown page is generated using:
 #src Literate.markdown("tutorials/Tutorial_AlpineData.jl","docs/src/man",keepcomments=true, execute=false, codefence = "```julia" => "```")
diff --git a/tutorials/Tutorial_Basic.jl b/tutorials/Tutorial_Basic.jl
index 8dfe9f90..6d4ec9d7 100644
--- a/tutorials/Tutorial_Basic.jl
+++ b/tutorials/Tutorial_Basic.jl
@@ -21,9 +21,17 @@ Topo_Alps = load_GMG("https://zenodo.org/records/10738510/files/AlpsTopo.jld2?do
 # We can write this to disk as well
 write_paraview(Topo_Alps,"Topo_Alps")
 
-# If we open both datasets in Paraview, we see this (after giving some color to the topography):
+# If we open both datasets in Paraview, and changing both files from outline/solid colors to the corresponding data field, we see:
+# ![Basic_Tutorial_1](../assets/img/Basic_Tutorial_Paraview_1.png)
+# Now we can change the colormap on the right side, marked by a red square. For topography we use the `Oleron` colormap, which you can download [here](https://www.fabiocrameri.ch/colourmaps/).
+# For the tomography we use the `Roma` scientific colormap. You will now see a blue'ish box of the tomography, this is not the best color to visualise the data. Let's invert the colormap by clicking on the item marked by the blue arrow.
+# Now we see the tomography in a more intuitive way, but the topography is not visible anymore. We can change the opacity of the tomography by setting a value in the `Opacity` field marked by the red square.
+# Note that you will need to adapt the range of the topography colormap as the change in color is not at 0.0. By clicking on the item marked by the black arrow, you can set your desired range.
+
+# ![Basic_Tutorial_1](../assets/img/Basic_Tutorial_Paraview_2.png)
+
+# Now you should see something like this:
 # ![Basic_Tutorial_1](../assets/img/Basic_Tutorial_1.png)
-# Note that I use the `Oleron` scientific colormap for the tomography which you can download [here](https://www.fabiocrameri.ch/colourmaps/)
 #
 # ### 2. Extract subset of data
 # As you can see the tomographic data covers a much larger area than the Alps itself, and in most of that area there is no data.
@@ -32,7 +40,7 @@ Tomo_Alps = extract_subvolume(Tomo_Alps_full,Lon_level=(4,20),Lat_level=(36,50),
 
 write_paraview(Tomo_Alps,"Tomo_Alps");
 
-# Which looks like:
+# After loading the new data again in paraview, switching to the proper data field and adjusting the colormap, you should see something like this:
 # ![Basic_Tutorial_2](../assets/img/Basic_Tutorial_2.png)
 
 # ### 3. Create cross sections
@@ -45,40 +53,40 @@ data_200km = cross_section(Tomo_Alps, Depth_level=-200)
 data_200km_exact = cross_section(Tomo_Alps, Depth_level=-200, Interpolate=true)
 
 # In general, you can get help info for all functions with `?`:
-# ```julia
+# ```julia-repl
 # help?> cross_section
 # search: cross_section cross_section_volume cross_section_points cross_section_surface flatten_cross_section
-# 
+#
 #   cross_section(DataSet::AbstractGeneralGrid; dims=(100,100), Interpolate=false, Depth_level=nothing, Lat_level=nothing, Lon_level=nothing, Start=nothing, End=nothing, Depth_extent=nothing, section_width=50km)
-# 
+#
 #   Creates a cross-section through a GeoData object.
-# 
+#
 #     •  Cross-sections can be horizontal (map view at a given depth), if Depth_level is specified
-# 
+#
 #     •  They can also be vertical, either by specifying Lon_level or Lat_level (for a fixed lon/lat), or by defining both Start=(lon,lat) & End=(lon,lat) points.
-# 
+#
 #     •  Depending on the type of input data (volume, surface or point data), cross sections will be created in a different manner:
-# 
+#
 #     1. Volume data: data will be interpolated or directly extracted from the data set.
-# 
+#
 #     2. Surface data: surface data will be interpolated or directly extracted from the data set
-# 
+#
 #     3. Point data: data will be projected to the chosen profile. Only data within a chosen distance (default is 50 km) will be used
-# 
+#
 #     •  Interpolate indicates whether we want to simply extract the data from the data set (default) or whether we want to linearly interpolate it on a new grid, which has dimensions as specified in dims NOTE: THIS ONLY APPLIES TO VOLUMETRIC AND SURFACE DATA
 #        SETS
-# 
+#
 #     •  'section_width' indicates the maximal distance within which point data will be projected to the profile
-# 
+#
 #   Example:
 #   ≡≡≡≡≡≡≡≡
-# 
+#
 #   julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
 #   julia> Data            =   Depth*2;                # some data
 #   julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
-#   julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))); 
-#   julia> Data_cross      =   cross_section(Data_set3D, Depth_level=-100km)  
-#   GeoData 
+#   julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)));
+#   julia> Data_cross      =   cross_section(Data_set3D, Depth_level=-100km)
+#   GeoData
 #     size  : (11, 11, 1)
 #     lon   ϵ [ 10.0 : 20.0]
 #     lat   ϵ [ 30.0 : 40.0]
@@ -93,8 +101,13 @@ Cross_vert = cross_section(Tomo_Alps, Start=(5,47), End=(15,44))
 write_paraview(Cross_vert,"Cross_vert");
 write_paraview(data_200km,"data_200km");
 
+# After loading the data in Paraview, you can use the `Clip` tool on the topography to only show the topography above sealevel and make it 60% transparent. Also adjust the colormap of the tomography to 5.0 and -5.0
+
+# ![Basic_Tutorial_3](../assets/img/Basic_Tutorial_Paraview_3.png)
+
+# After doing all these steps, you should see something like this:
 # ![Basic_Tutorial_3](../assets/img/Basic_Tutorial_3.png)
-# In creating this image, I used the `Clip` tool of Paraview to only show topography above sealevel and made it 50% transparent.
+
 #
 # ### 4. Cartesian data
 # As you can see, the curvature or the Earth is taken into account here. Yet, for many applications it is more convenient to work in Cartesian coordinates (kilometers) rather then in geographic coordinates.
@@ -120,13 +133,13 @@ Yet, because of the curvature of the Earth, the resulting 3D model is not strict
 This can be achieved in a relatively straightforward manner, by creating a new 3D dataset that is slightly within the curved boundaries of the projected data set:
 =#
 
-Tomo_rect = CartData(xyz_grid(-550.0:10:600, -500.0:10:700, -600.0:5:-17)); 
+Tomo_rect = CartData(xyz_grid(-550.0:10:600, -500.0:10:700, -600.0:5:-17));
 #
 # the routine `project_CartData` will then project the data from the geographic coordinates to the new rectilinear grid:
 Tomo_rect = project_CartData(Tomo_rect, Tomo_Alps, proj)
 
 # we can do the same with topography:
-Topo_rect = CartData(xyz_grid(-550.0:1:600, -500.0:1:700, 0)) 
+Topo_rect = CartData(xyz_grid(-550.0:1:600, -500.0:1:700, 0))
 Topo_rect = project_CartData(Topo_rect, Topo_Alps, proj)
 
 # Save it:
@@ -134,9 +147,8 @@ write_paraview(Tomo_rect,"Tomo_rect");
 write_paraview(Topo_rect,"Topo_rect");
 # ![Basic_Tutorial_5](../assets/img/Basic_Tutorial_5.png)
 
-# At this stage, the data can directly be used to generate cartesian numerical model setups, as explained in the other tutorials. 
+# At this stage, the data can directly be used to generate cartesian numerical model setups, as explained in the other tutorials.
 
 
 #src Note: The markdown page is generated using:
 #src Literate.markdown("tutorials/Tutorial_Basic.jl","docs/src/man",keepcomments=true, execute=true, codefence = "```julia" => "```")
-
diff --git a/tutorials/Tutorial_Chmy_MPI.jl b/tutorials/Tutorial_Chmy_MPI.jl
index 6e6de8b2..f1cb2178 100644
--- a/tutorials/Tutorial_Chmy_MPI.jl
+++ b/tutorials/Tutorial_Chmy_MPI.jl
@@ -3,8 +3,8 @@
 # ## Aim
 # In this tutorial, your will learn how to use [Chmy](https://github.com/PTsolvers/Chmy.jl) to perform a 2D diffusion simulation
 # on one or multiple CPU's or GPU's.
-# `Chmy` is a package that allows you to specify grids and fields and create finite difference simulations, by 
-# applying stencil-based kernels in an efficient manner.   
+# `Chmy` is a package that allows you to specify grids and fields and create finite difference simulations, by
+# applying stencil-based kernels in an efficient manner.
 #
 # ## 1. Load Chmy and required packages
 using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.Fields, Chmy.BoundaryConditions, Chmy.GridOperators, Chmy.KernelLaunch
@@ -54,7 +54,7 @@ end
     dims_g = dims_l .* dims(topo)
     grid   = UniformGrid(arch; origin=(-2, -2), extent=(4, 4), dims=dims_g)
     launch = Launcher(arch, grid, outer_width=(16, 8))
-    
+
     ##@info "mpi" me grid
 
     nx, ny = dims_g
@@ -65,17 +65,17 @@ end
     ## allocate fields
     C = Field(backend, grid, Center())
     P = Field(backend, grid, Center(), Int32)   # phases
-    
+
     q = VectorField(backend, grid)
     C_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(C)) .* dims(topo)) : nothing
-    
-    ## Use the `GeophysicalModelGenerator` to set the initial conditions. Note that 
+
+    ## Use the `GeophysicalModelGenerator` to set the initial conditions. Note that
     ## you have to call this for a `Phases` and a `Temp` grid, which we call `C` here.
     add_box!(P,C,grid,  xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400))
 
     ## set BC's and updates the halo:
     bc!(arch, grid, C => Neumann(); exchange=C)
-    
+
     ## visualisation
     fig = Figure(; size=(400, 320))
     ax  = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y", title="it = 0")
@@ -102,7 +102,7 @@ end
 
 # In the code above, the part that calls `GMG` is:
 
-# ```julia 
+# ```julia
 # add_box!(P,C,grid,  xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400))
 # ```
 # which works just like any of the other GMG function, except that you pass Chmy Scalar Fields and a Chmy grid object.
@@ -131,7 +131,7 @@ MPI.Finalize()
 # mpiexecjl -n 4 --project=. julia Tutorial_Chmy_MPI.jl
 # ```
 
-# The full file can be downloaded [here](../../../tutorials/Tutorial_Chmy_MPI.jl)
+# The full file can be downloaded [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorials/Tutorial_Chmy_MPI.jl)
 
 #src Note: The markdown page is generated using:
 #src Literate.markdown("tutorials/Tutorial_Chmy_MPI.jl","docs/src/man",keepcomments=true, execute=false, codefence = "```julia" => "```")