diff --git a/.github/workflows/blank.yml b/.github/workflows/CI.yml
similarity index 100%
rename from .github/workflows/blank.yml
rename to .github/workflows/CI.yml
diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
new file mode 100644
index 00000000..da6ab6df
--- /dev/null
+++ b/.github/workflows/Documenter.yml
@@ -0,0 +1,48 @@
+# Note: this workflow is copied from the Trixi.jl package
+name: Documentation
+on:
+  push:
+    branches:
+      - 'main'
+    tags: '*'
+    paths-ignore:
+      - '.zenodo.json'
+      - '.github/workflows/benchmark.yml'
+      - '.github/workflows/CI.yml'
+      - '.github/workflows/CompatHelper.yml'
+  pull_request:
+    paths-ignore:
+      - '.zenodo.json'
+      - '.github/workflows/CI.yml'
+      - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
+      - 'benchmark/**'
+      - 'utils/**'
+  workflow_dispatch:
+
+# Cancel redundant CI tests automatically
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1.9'
+          show-versioninfo: true
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON: ""
+      - name: Install dependencies
+        run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+          GKSwstype: 100 # To make GitHub Action work, disable showing a plot window with the GR backend of the Plots package
+        run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs --color=yes docs/make.jl
\ No newline at end of file
diff --git a/.github/workflows/cleanup_caches.yml b/.github/workflows/cleanup_caches.yml
new file mode 100644
index 00000000..8985f34c
--- /dev/null
+++ b/.github/workflows/cleanup_caches.yml
@@ -0,0 +1,29 @@
+name: cleanup caches by a branch
+on:
+  pull_request:
+    types:
+      - closed
+
+jobs:
+  cleanup:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Cleanup
+        run: |
+          gh extension install actions/gh-actions-cache
+
+          echo "Fetching list of cache key"
+          cacheKeysForPR=$(gh actions-cache list -R $REPO -B $BRANCH -L 100 | cut -f 1 )
+
+          ## Setting this to not fail the workflow while deleting cache keys.
+          set +e
+          echo "Deleting caches..."
+          for cacheKey in $cacheKeysForPR
+          do
+              gh actions-cache delete $cacheKey -R $REPO -B $BRANCH --confirm
+          done
+          echo "Done"
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          REPO: ${{ github.repository }}
+          BRANCH: refs/pull/${{ github.event.pull_request.number }}/merge
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
deleted file mode 100644
index 0312a351..00000000
--- a/.github/workflows/documentation.yml
+++ /dev/null
@@ -1,24 +0,0 @@
-name: Documentation
-
-on:
-  push:
-    branches:
-      - main
-    tags: '*'
-  pull_request:
-
-jobs:
-  build:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@latest
-        with:
-          version: '1.9'
-      - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.build(); Pkg.develop(PackageSpec(path=pwd()));  Pkg.instantiate()'
-      - name: Build and deploy
-        env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
diff --git a/docs/.gitignore b/docs/.gitignore
new file mode 100644
index 00000000..aa3e5d8c
--- /dev/null
+++ b/docs/.gitignore
@@ -0,0 +1 @@
+gmt.history
\ No newline at end of file
diff --git a/docs/Project.toml b/docs/Project.toml
index a841a404..cdd27916 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,8 @@
 [deps]
-GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
+GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 
 [compat]
-Documenter = "0.26"
+Documenter = "1"
diff --git a/docs/gmt.history b/docs/gmt.history
index b53e2c35..2d903eea 100644
--- a/docs/gmt.history
+++ b/docs/gmt.history
@@ -3,6 +3,6 @@ BEGIN GMT 6.5.0
 B	afWSen
 J	X
 JX	X15c/0
-R	0/1/0.1/1
+R	0/1/0/0.9
 @L	1
 END
diff --git a/docs/make.jl b/docs/make.jl
index 5174d41f..828e253e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,7 +1,12 @@
-using GeophysicalModelGenerator
-using .GeophysicalModelGenerator
 using Documenter
 
+push!(LOAD_PATH, dirname(@__DIR__))
+
+using GeophysicalModelGenerator
+
+# Importing these activates package extensions
+import GLMakie, GMT
+
 #DocMeta.setdocmeta!(GeophysicalModelGenerator, :DocTestSetup, :(using GeophysicalModelGenerator); recursive=true)
 
 # Get GeophysicalModelGenerator.jl root directory
@@ -48,6 +53,7 @@ open(joinpath(@__DIR__, "src", "man", "code_of_conduct.md"), "w") do io
     println(io, "> ", line)
   end
 end
+
 open(joinpath(@__DIR__, "src", "man", "contributing.md"), "w") do io
     # Point to source license file
     println(io, """
@@ -67,7 +73,7 @@ open(joinpath(@__DIR__, "src", "man", "contributing.md"), "w") do io
 
 makedocs(;
     modules=[GeophysicalModelGenerator],
-    authors="Marcel Thielmann, Boris Kaus",
+    authors="Boris Kaus, Marcel Thielmann",
     sitename="GeophysicalModelGenerator.jl",
     format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
@@ -76,21 +82,21 @@ makedocs(;
         "Home" => "index.md",
         "Tutorials" => Any[
             "Overview" =>  "man/tutorials.md",
-            "3D seismic tomography from ASCII" =>  "man/tutorial_load3DSeismicData.md",
-            "3D seismic tomography from netCDF" =>  "man/tutorial_loadregular3DSeismicData_netCDF.md",
-            "Visualize Moho topography" =>  "man/tutorial_MohoTopo.md",
-            "Create GMT-based topography" =>  "man/tutorial_GMT_Topography.md",
-            "Coastlines" =>  "man/tutorial_Coastlines.md",
-            "Import screenshots" =>  "man/tutorial_Screenshot_To_Paraview.md",
-            "Interpolate irregular 3D seismic tomography" =>  "man/tutorial_loadirregular3DSeismicData.md",
-            "ETOPO1 Topography and geological maps" =>  "man/tutorial_GMT_Topography_GeologicalMap.md",
-            "ISC earthquake data" =>  "man/tutorial_ISC_data.md",
-            "Plot GPS vectors" =>  "man/tutorial_GPS.md",
-            "Read UTM data" =>  "man/tutorial_UTM.md",
-            "VoteMaps" =>  "man/Tutorial_Votemaps.md",
-            "Kilometer-scale volcano" =>  "man/tutorial_local_Flegrei.md",
-            "Generating LaMEM model" =>  "man/LaPalma_example.md",
-            "Create movies" =>  "man/tutorial_time_Seismicity.md"
+            "1 - 3D seismic tomography from ASCII" =>  "man/tutorial_load3DSeismicData.md",
+            "2 - 3D seismic tomography from netCDF" =>  "man/tutorial_loadregular3DSeismicData_netCDF.md",
+            "3 - Visualize Moho topography" =>  "man/tutorial_MohoTopo.md",
+            "4 - Create GMT-based topography" =>  "man/tutorial_GMT_Topography.md",
+            "5 - Coastlines" =>  "man/tutorial_Coastlines.md",
+            "6 - Import screenshots" =>  "man/tutorial_Screenshot_To_Paraview.md",
+            "7 - Interpolate irregular 3D seismic tomography" =>  "man/tutorial_loadirregular3DSeismicData.md",
+            "8 - ETOPO1 Topography and geological maps" =>  "man/tutorial_GMT_Topography_GeologicalMap.md",
+            "9 - ISC earthquake data" =>  "man/tutorial_ISC_data.md",
+            "10 - Plot GPS vectors" =>  "man/tutorial_GPS.md",
+            "11 - Read UTM data" =>  "man/tutorial_UTM.md",
+            "12 - VoteMaps" =>  "man/Tutorial_Votemaps.md",
+            "13 - Campi Flegrei" =>  "man/tutorial_local_Flegrei.md",
+            "14 - Cartesian Volcano Model" =>  "man/LaPalma_example.md",
+            "15 - Create movies" =>  "man/tutorial_time_Seismicity.md"
         ],
         "User Guide" => Any[
             "Installation" =>  "man/installation.md",
@@ -110,6 +116,8 @@ makedocs(;
         "Code of Conduct" => "man/code_of_conduct.md",
         "License" => "man/license.md"
     ],
+    pagesonly=true,
+    warnonly=true
 )
 
 deploydocs(;
diff --git a/docs/src/man/LaPalma_example.md b/docs/src/man/LaPalma_example.md
index aaf34c79..cfa201d0 100644
--- a/docs/src/man/LaPalma_example.md
+++ b/docs/src/man/LaPalma_example.md
@@ -1,11 +1,12 @@
 ```@meta
-EditURL = "<unknown>/LaPalma_example.jl"
+EditURL = "../../../tutorials/LaPalma_example.jl"
 ```
 
-# Create a 3D LaMEM setup for La Palma
+# 14 - Create a Cartesian Model Setup for La Palma
 
 ## Aim
-In this tutorial, your will learn how to use real data to create a geodynamic model setup with LaMEM. We will use the data of La Palma, which is a volcanic island that started erupting in mid september 2021.
+In this tutorial, your will learn how to use real data to create a geodynamic model setup with LaMEM.
+We will use the data of the Cumbre Viejo eruption in La Palma, which is a volcanic island that erupted from mid september 2021 - december 2021.
 LaMEM is a cartesian geodynamic model, which implies that you will have to transfer the data from `GeoData` to `CartData`.
 
 
@@ -14,69 +15,44 @@ We will use two types of data to create the model
  1) Topography
  2) Earthquake locations
 
-We start with loading the required packages
+We start with loading the required packages, which includes `GMT` to download topography (an optional dependency for `GeophysicalModelGenerator`)
 
-````julia
-using GeophysicalModelGenerator, JLD2
-using GMT
-````
+```julia
+using GeophysicalModelGenerator, GMT
+```
 
 In case you managed to install GMT on your machine, you can automatically download the topography with
 
-````julia
+```julia
 Topo = ImportTopo(lon = [-18.7, -17.1], lat=[28.0, 29.2], file="@earth_relief_03s.grd")
-````
-
-````
-GeoData
-  size  : (1921, 1441, 1)
-  lon   ϵ [ 341.3 : 342.9]
-  lat   ϵ [ 28.0 : 29.2]
-  depth ϵ [ -4.38297802734375 km : 2.414 km]
-  fields: (:Topography,)
-````
-
-!!! note
-
-    In case you did not manage that, we have prepared a JLD2 file [here](https://seafile.rlp.net/f/64e0b0a6b1e941b7bba5/?dl=1). Download it, and doublecheck that you are in the same directory as the data file with:
-    ````julia
-      pwd()
-      "/Users/kausb/.julia/dev/GeophysicalModelGenerator/docs/src/scripts"
-    ````
-    Load the data with:
-    ````julia
-      julia> Topo=load_GMG("Topo_LaPalma")
-      GeoData
-        size      : (1921, 1441, 1)
-        lon       ϵ [ 341.3 : 342.9]
-        lat       ϵ [ 28.0 : 29.2]
-        depth     ϵ [ -4.090296875 : 0.0]
-        fields    : (:Topography,)
-        attributes: ["note"]
-    ````
+```
+
+In case you did not manage that, we have prepared a JLD2 file [here](https://seafile.rlp.net/f/03286a46a9f040a69b0f/?dl=1).
+Download it, and doublecheck that you are in the same directory as the data file with:
+
+```julia
+pwd()
+```
+
+Load the data:
+
+```julia
+Topo=load("Topo_LaPalma.jld2","Topo")
+```
 
 Next, lets load the seismicity. We have prepared a file with earthquake locations up to early November (from [https://www.ign.es/web/ign/portal/vlc-catalogo](https://www.ign.es/web/ign/portal/vlc-catalogo)).
-Download that [here](https://seafile.rlp.net/f/64e0b0a6b1e941b7bba5/?dl=1)
-
-````julia
-data_all_EQ = load_GMG("EQ_Data")
-GeoData
-  size      : (6045,)
-  lon       ϵ [ -17.841 : -17.8268]
-  lat       ϵ [ 28.5717 : 28.573]
-  depth     ϵ [ -29.0 : -11.0]
-  fields    : (:Magnitude, :TimeSinceJan1_2021)
-  attributes: ["note"]
-````
+Download that [here](https://seafile.rlp.net/f/03286a46a9f040a69b0f/?dl=1)
+
+```julia
+data_all_EQ = load("EQ_Data.jld2","data_all_EQ")
+```
 
 Write the data to paraview with
-````julia
-julia> Write_Paraview(data_all_EQ,"data_all_EQ",PointsData=true)
-Saved file: data_all_EQ.vtu
 
-julia> Write_Paraview(Topo,"Topo")
-Saved file: Topo.vts
-````
+```julia
+Write_Paraview(data_all_EQ,"data_all_EQ",PointsData=true)
+Write_Paraview(Topo,"Topo")
+```
 
 As earthquakes are point-wise data, you have to specify this.
 ![LaPalma_EQTopo_GeoData](../assets/img/TopoEQs_LaPalma_GeoData.png)
@@ -86,55 +62,41 @@ Note that this data is not in "easy" coordinates (see coordinate axis in the plo
 In order to create model setups, it is helpful to first transfer the data to Cartesian.
 This requires us to first determine a *projection point*, that is fixed. Often, it is helpful to use the center of the topography for this. In the present example, we will center the model around La Palma itself:
 
-````julia
-julia> proj = ProjectionPoint(Lon=-17.84, Lat=28.56)
-ProjectionPoint(28.56, -17.84, 222158.69478000276, 3.1625327286383654e6, 28, true)
-````
+```julia
+proj = ProjectionPoint(Lon=-17.84, Lat=28.56)
+```
 
 Once this is done you can convert the topographic data to the cartesian reference frame
 
-````julia
-julia> EQ_cart   = Convert2CartData(data_all_EQ, proj);
-julia> Topo_cart = Convert2CartData(Topo, proj)
-CartData
-    size    : (1921, 1441, 1)
-    x       ϵ [ -86.09445705828863 : 73.67229892155609]
-    y       ϵ [ -63.5531883197492 : 73.28446155584604]
-    z       ϵ [ -4.38352294921875 : 2.414]
-    fields  : (:Topography,)
-  attributes: ["note"]
-````
+```julia
+EQ_cart   = Convert2CartData(data_all_EQ, proj);
+Topo_cart = Convert2CartData(Topo, proj)
+```
 
 It is important to realize that the cartesian coordinates of the topographic grid is no longer strictly orthogonal after this conversion. You don't notice that in the current example, as the model domain is rather small.
 In other cases, however, this is quite substantial (e.g., India-Asia collision zone).
 LaMEM needs an orthogonal grid of topography, which we can create with:
 
-````julia
-julia> Topo_LaMEM = CartData(XYZGrid(-70:.2:70,-60:.2:70,0));
-````
+```julia
+Topo_LaMEM = CartData(XYZGrid(-70:.2:70,-60:.2:70,0));
+nothing #hide
+```
 
 In a next step, the routine `ProjectCartData` projects a `GeoData` structure to a `CartData` struct
 
-````julia
-julia> Topo_LaMEM = ProjectCartData(Topo_LaMEM, Topo, proj)
-CartData
-    size    : (701, 651, 1)
-    x       ϵ [ -70.0 : 70.0]
-    y       ϵ [ -60.0 : 70.0]
-    z       ϵ [ -4.29608214262314 : 2.3840784607152257]
-    fields  : (:Topography,)
-  attributes: ["note"]
-````
+```julia
+Topo_LaMEM = ProjectCartData(Topo_LaMEM, Topo, proj)
+```
 
 Let's have a look at the data:
 
-````julia
-julia> Write_Paraview(EQ_cart,"EQ_cart",PointsData=true)
-Saved file: EQ_cart.vtu
+```julia
+Write_Paraview(EQ_cart,"EQ_cart",PointsData=true)
+Write_Paraview(Topo_LaMEM,"Topo_LaMEM")
 
-julia>  Write_Paraview(Topo_LaMEM,"Topo_LaMEM")
-Saved file: Topo_LaMEM.vts
-````
+
+#=
+```
 
 ## 3. Create LaMEM setup
 
@@ -143,62 +105,61 @@ the number of control volumes (elements), you want to apply in every direction.
 The LaMEM input file can be downloaded [here](../scripts/LaPalma.dat).
 Make sure you are in the same directory as the *.dat file & execute the following command
 
-````julia
-julia> Grid    = ReadLaMEM_InputFile("LaPalma.dat")
-LaMEM Grid:
-  nel         : (64, 64, 32)
-  marker/cell : (3, 3, 3)
-  markers     : (192, 192, 96)
-  x           ϵ [-50.0 : 50.0]
-  y           ϵ [-40.0 : 40.0]
-  z           ϵ [-80.0 : 15.0]
-````
+```julia
+Grid    = ReadLaMEM_InputFile("LaPalma.dat")
+```
 
 The `LaMEM_grid` structure contains the number of elements in every direction and the number of markers in every cell.
 It also contains `Grid.X`, `Grid.Y`, `Grid.Z`, which are the coordinates of each of the markers in the 3 directions.
 
 In a next step we need to give each of these points a `Phase` number (which is an integer, that indicates the type of the rock that point has), as well as the temperature (in Celsius).
 
-````julia
-julia> Phases  = ones(Int32,size(Grid.X))*2;
-julia> Temp    = ones(Float64,size(Grid.X));
-````
+```julia
+Phases  = ones(Int32,size(Grid.X))*2;
+Temp    = ones(Float64,size(Grid.X));
+nothing #hide
+```
 
 In this example we set the temperature based on the depth, assuming a constant geotherm. Depth is given in kilometers
 
-````julia
-julia> Geotherm =  30;
-julia> Temp     =  -Grid.Z.*Geotherm;
-````
+```julia
+Geotherm =  30;
+Temp     =  -Grid.Z.*Geotherm;
+nothing #hide
+```
 
 Since we are in La Palma, we assume that temperatures are not less than 20 degrees. Moreover, in the mantle, we have the mantle adiabat (1350 C)
 
-````julia
-julia> Temp[Temp.<20]    .=  20;
-julia> Temp[Temp.>1350]  .=  1350;
-````
+```julia
+Temp[Temp.<20]    .=  20;
+Temp[Temp.>1350]  .=  1350;
+nothing #hide
+```
 
 Because of lacking data, we have pretty much no idea where the Moho is in La Palma.
 Somewhat arbitrary, we assume it to be at 40 km depth, and set the rocks below that to "mantle":
 
-````julia
-julia> ind = findall( Grid.Z .< -40);
-julia> Phases[ind] .= 3;
-````
+```julia
+ind = findall( Grid.Z .< -40);
+Phases[ind] .= 3;
+nothing #hide
+```
 
 Everything above the free surface is assumed to be "air"
 
-````julia
-julia> ind     =   AboveSurface(Grid, Topo_cart);
-julia> Phases[ind] .= 0;
-````
+```julia
+ind     =   AboveSurface(Grid, Topo_cart);
+Phases[ind] .= 0;
+nothing #hide
+```
 
 And all "air" points that are below sea-level becomes "water"
 
-````julia
-julia> ind = findall( (Phases.==0) .& (Grid.Z .< 0));
-julia> Phases[ind] .= 1;
-````
+```julia
+ind = findall( (Phases.==0) .& (Grid.Z .< 0));
+Phases[ind] .= 1;
+nothing #hide
+```
 
 You can interpret the seismicity in different ways. If you assume that there are fully molten magma chambers, there should be no earthquakes within the magma chamber (as stresses are liklely to be very small).
 If, however, this is a partially molten mush, extraction of melt of that mush will change the fluid pressure and may thus release build-up stresses.
@@ -206,46 +167,50 @@ We will assume the second option in our model setup.
 
 Looking at the seismicity, there is a swarm at around 35 km depth
 
-````julia
-julia> AddSphere!(Phases,Temp,Grid, cen=(0,0,-35), radius=5,phase=ConstantPhase(5), T=ConstantTemp(1200));
-````
+```julia
+AddSphere!(Phases,Temp,Grid, cen=(0,0,-35), radius=5, phase=ConstantPhase(5), T=ConstantTemp(1200));
+nothing #hide
+```
 
 A shallower one exists as well
 
-````julia
-julia> AddEllipsoid!(Phases,Temp,Grid, cen=(-1,0,-11), axes=(3,3,8), StrikeAngle=225, DipAngle=45, phase=ConstantPhase(5), T=ConstantTemp(1200));
-````
+```julia
+AddEllipsoid!(Phases,Temp,Grid, cen=(-1,0,-11), axes=(3,3,8), StrikeAngle=225, DipAngle=45, phase=ConstantPhase(5), T=ConstantTemp(1200));
+nothing #hide
+```
 
 And there might be a mid-crustal one
 
-````julia
-julia> AddEllipsoid!(Phases,Temp,Grid, cen=(-0,0,-23), axes=(8,8,2), StrikeAngle=0, DipAngle=0, phase=ConstantPhase(5), T=ConstantTemp(1200));
-````
+```julia
+AddEllipsoid!(Phases,Temp,Grid, cen=(-0,0,-23), axes=(8,8,2), StrikeAngle=0, DipAngle=0, phase=ConstantPhase(5), T=ConstantTemp(1200));
+nothing #hide
+```
 
 We can generate a 3D model setup, which must include the `Phases` and `Temp` arrays
 
-````julia
-julia> Model3D = CartData(Grid, (Phases=Phases,Temp=Temp));
-julia> Write_Paraview(Model3D,"Model3D")
-Saved file: Model3D.vts
-````
+```julia
+Model3D = CartData(Grid, (Phases=Phases,Temp=Temp))
+Write_Paraview(Model3D,"Model3D")
+```
 
 The model setup looks like this
 ![LaMEM_ModelSetup_LaPalma](../assets/img/LaMEM_ModelSetup_LaPalma.png)
 
 We can create a LaMEM marker file from the `Model3D` setup and the (cartesian) topography
 
-````julia
-julia> Save_LaMEMTopography(Topo_cart, "Topography.txt")
-Written LaMEM topography file: Topography.txt
-julia> Save_LaMEMMarkersParallel(Model3D)
-Writing LaMEM marker file -> ./markers/mdb.00000000.dat
-````
-
+```julia
+Save_LaMEMTopography(Topo_cart, "Topography.txt")
+Save_LaMEMMarkersParallel(Model3D)
+```
 
 Next, you can use this to run the LaMEM model and visualize the model results in Paraview as well.
 If you are interested in doing this, have a look at the LaMEM [wiki](https://bitbucket.org/bkaus/lamem/wiki/Home) pages.
 
+```julia
+=#
+```
+
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+
diff --git a/docs/src/man/Tutorial_Votemaps.md b/docs/src/man/Tutorial_Votemaps.md
index 2d7e4078..ce05f3b1 100644
--- a/docs/src/man/Tutorial_Votemaps.md
+++ b/docs/src/man/Tutorial_Votemaps.md
@@ -1,8 +1,8 @@
 ```@meta
-EditURL = "<unknown>/Tutorial_Votemaps.jl"
+EditURL = "../../../tutorials/Tutorial_Votemaps.jl"
 ```
 
-# Votemaps
+# 12 - Votemaps
 
 ## Aim
 In this tutorial, your will learn how to create Votemaps that compare different tomographic models and look for similarities between different models.
@@ -58,10 +58,10 @@ The result of this gives a feeling which features are consistent between the 3 m
 So how do we create Votemaps?
 Doing this is rather simple:
 
-````@example Tutorial_Votemaps
+```julia
 Data_VoteMap = VoteMap( [Pwave_Paffrath,       PSwave_Koulakov,    Pwave_Zhao],
                         ["dVp_Percentage>3.0","dVp_percentage>2.0","dVp_Percentage>2.0"], dims=(100,100,100))
-````
+```
 
 This will look at the common `lon`,`lat`,`depth` ranges between all 3 models, interpret each of the models to a common grid of size `(100,100,100)` and apply each of the criteria specified
 The resulting `GeoData` struct looks like:
diff --git a/docs/src/man/datastructures.md b/docs/src/man/datastructures.md
index 8f1e8983..3843408d 100644
--- a/docs/src/man/datastructures.md
+++ b/docs/src/man/datastructures.md
@@ -5,11 +5,13 @@ We also provide a `UTMData`, which is essentially the same but with UTM coordina
 For plotting, we transfer this into the `ParaviewData` structure, which has cartesian coordinates around the center of the Earth. We employ the `wgs84` reference ellipsoid as provided by the [Geodesy.jl](https://github.com/JuliaGeo/Geodesy.jl) package to perform this transformation.
 
 ```@docs
-GeophysicalModelGenerator.GeoData
-GeophysicalModelGenerator.UTMData
-GeophysicalModelGenerator.ParaviewData
-GeophysicalModelGenerator.CartData
-GeophysicalModelGenerator.LonLatDepthGrid
-GeophysicalModelGenerator.XYZGrid
-GeophysicalModelGenerator.ProjectionPoint
+UTMData
+ParaviewData
+CartData
+LonLatDepthGrid
+XYZGrid
+ProjectionPoint
 ```
+
+
+GeoData
diff --git a/docs/src/man/lamem.md b/docs/src/man/lamem.md
index 4f63e9a2..37d89251 100644
--- a/docs/src/man/lamem.md
+++ b/docs/src/man/lamem.md
@@ -25,6 +25,7 @@ GeophysicalModelGenerator.ConstantTemp
 GeophysicalModelGenerator.LinearTemp
 GeophysicalModelGenerator.HalfspaceCoolingTemp
 GeophysicalModelGenerator.SpreadingRateTemp
+GeophysicalModelGenerator.LithosphericTemp
 GeophysicalModelGenerator.ConstantPhase
 GeophysicalModelGenerator.Compute_Phase
 GeophysicalModelGenerator.LithosphericPhases
diff --git a/docs/src/man/listfunctions.md b/docs/src/man/listfunctions.md
index 022c6c2e..c1b3b8d6 100644
--- a/docs/src/man/listfunctions.md
+++ b/docs/src/man/listfunctions.md
@@ -1,6 +1,7 @@
 # List of all functions
 
 Here an overview of all functions:
-```@index
+```@autodocs
+Modules = [GeophysicalModelGenerator]
 ```
 
diff --git a/docs/src/man/projection.md b/docs/src/man/projection.md
index eee26316..4b777735 100644
--- a/docs/src/man/projection.md
+++ b/docs/src/man/projection.md
@@ -102,8 +102,10 @@ So this interpolates the topographic data from the `GeoData` to the orthogonal c
 You can do similar projections with full 3D data sets or pointwise data. 
 
 #### 3. List of functions
-```@docs
-GeophysicalModelGenerator.Convert2CartData
-GeophysicalModelGenerator.ProjectCartData
-GeophysicalModelGenerator.Convert2UTMzone
-```
+
+
+# ```@docs
+# GeophysicalModelGenerator.Convert2CartData
+# GeophysicalModelGenerator.ProjectCartData
+# GeophysicalModelGenerator.Convert2UTMzone
+# ```
diff --git a/docs/src/man/tools.md b/docs/src/man/tools.md
index 2d17c007..7b5d2361 100644
--- a/docs/src/man/tools.md
+++ b/docs/src/man/tools.md
@@ -11,6 +11,7 @@ GeophysicalModelGenerator.SubtractHorizontalMean
 GeophysicalModelGenerator.AboveSurface
 GeophysicalModelGenerator.BelowSurface
 GeophysicalModelGenerator.InterpolateDataOnSurface
+GeophysicalModelGenerator.InterpolateTopographyOnPlane
 GeophysicalModelGenerator.ParseColumns_CSV_File
 GeophysicalModelGenerator.RotateTranslateScale!
 GeophysicalModelGenerator.Convert2UTMzone
diff --git a/docs/src/man/tutorial_Coastlines.md b/docs/src/man/tutorial_Coastlines.md
index 39c39ad2..675e9551 100644
--- a/docs/src/man/tutorial_Coastlines.md
+++ b/docs/src/man/tutorial_Coastlines.md
@@ -1,4 +1,4 @@
-# Add coastlines 
+# 5 - Add coastlines 
 
 ## Goal
 
diff --git a/docs/src/man/tutorial_GMT_Topography.md b/docs/src/man/tutorial_GMT_Topography.md
index d233fb38..7e0514ce 100644
--- a/docs/src/man/tutorial_GMT_Topography.md
+++ b/docs/src/man/tutorial_GMT_Topography.md
@@ -1,4 +1,4 @@
-# Extract topographic data from GMT.jl 
+# 4 - Extract topographic data from GMT.jl 
 
 ## Goal
 
diff --git a/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md b/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md
index c8abcaa1..48535876 100644
--- a/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md
+++ b/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md
@@ -1,4 +1,4 @@
-# Extract ETOPO1 topographic data using GMT.jl and drape a geological map on top of the topography (given as raster graphics)
+# 8 - Extract ETOPO1 topographic data using GMT.jl and drape a geological map on top of the topography (given as raster graphics)
 
 ## Goal
 
diff --git a/docs/src/man/tutorial_ISC_data.md b/docs/src/man/tutorial_ISC_data.md
index 8d43dad9..07f6e9c6 100644
--- a/docs/src/man/tutorial_ISC_data.md
+++ b/docs/src/man/tutorial_ISC_data.md
@@ -1,4 +1,4 @@
-# Plot ISC earthquake data
+# 10 - Plot ISC earthquake data
 
 ## Goal
 This explains how to load earthquake data obtained from the ISC catalogue.
diff --git a/docs/src/man/tutorial_MohoTopo.md b/docs/src/man/tutorial_MohoTopo.md
index 64040ae6..16b607b2 100644
--- a/docs/src/man/tutorial_MohoTopo.md
+++ b/docs/src/man/tutorial_MohoTopo.md
@@ -1,4 +1,4 @@
-# Moho topography
+# 3 - Moho topography
 
 ## Goal
 This explains how to load the Moho topography for Italy and the Alps and create a paraview file
diff --git a/docs/src/man/tutorial_Screenshot_To_Paraview.md b/docs/src/man/tutorial_Screenshot_To_Paraview.md
index 7278c2a6..5375dcd9 100644
--- a/docs/src/man/tutorial_Screenshot_To_Paraview.md
+++ b/docs/src/man/tutorial_Screenshot_To_Paraview.md
@@ -1,4 +1,4 @@
-# Import profiles/maps from published papers
+# 6 - Import profiles/maps from published papers
 
 ## Goal
 Ideally, all data should be available in digital format, after which you could use the tools described in the other tutorial to transform them into `GeoData` and export them to VTK.
@@ -7,14 +7,14 @@ Yet, the reality is different and often data is not (yet) available, or papers a
 For that reason, `GeophysicalModelGenerator` has tools that allow you to transfer a screenshot from any published paper into `GeoData/Paraview` and see it in 3D at the correct geographic location. This can be done for vertical profiles and for mapviews, which gives you a quick and easy way to see those papers in a new (3D) light.
 
 Here, we explain how.
-- [Import profiles/maps from published papers](#import-profilesmaps-from-published-papers)
+- [6 - Import profiles/maps from published papers](#6---import-profilesmaps-from-published-papers)
   - [Goal](#goal)
   - [General procedure](#general-procedure)
       - [1. Download data and crop images](#1-download-data-and-crop-images)
-      - [2. Read data of a cross-section & create VTS file](#2-read-data-of-a-cross-section--create-vts-file)
-      - [3. Read data of a mapview & create *.vts file](#3-read-data-of-a-mapview--create-vts-file)
+      - [2. Read data of a cross-section \& create VTS file](#2-read-data-of-a-cross-section--create-vts-file)
+      - [3. Read data of a mapview \& create \*.vts file](#3-read-data-of-a-mapview--create-vts-file)
       - [4. Using an automatic digitizer to pick points on map](#4-using-an-automatic-digitizer-to-pick-points-on-map)
-      - [5. Creating a multiblock Paraview/*.vtm file](#5-creating-a-multiblock-paraviewvtm-file)
+      - [5. Creating a multiblock Paraview/\*.vtm file](#5-creating-a-multiblock-paraviewvtm-file)
       - [6. Julia script](#6-julia-script)
 ## General procedure
 #### 1. Download data and crop images
diff --git a/docs/src/man/tutorial_UTM.md b/docs/src/man/tutorial_UTM.md
index e13ef20d..f6fd1753 100644
--- a/docs/src/man/tutorial_UTM.md
+++ b/docs/src/man/tutorial_UTM.md
@@ -1,4 +1,4 @@
-# Read in UTM data
+# 11 - Read in UTM data
 
 ## Goal
 
diff --git a/docs/src/man/tutorial_load3DSeismicData.md b/docs/src/man/tutorial_load3DSeismicData.md
index 24909c1a..92d9991b 100644
--- a/docs/src/man/tutorial_load3DSeismicData.md
+++ b/docs/src/man/tutorial_load3DSeismicData.md
@@ -1,4 +1,4 @@
-# 3D tomography model in CSV formation
+# 1 - 3D tomography model in CSV formation
 
 ## 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;
diff --git a/docs/src/man/tutorial_loadirregular3DSeismicData.md b/docs/src/man/tutorial_loadirregular3DSeismicData.md
index 3a5575e1..fd0d21ec 100644
--- a/docs/src/man/tutorial_loadirregular3DSeismicData.md
+++ b/docs/src/man/tutorial_loadirregular3DSeismicData.md
@@ -1,4 +1,4 @@
-# 3D tomography model that is re-interpolated on a regular grid
+# 7 - 3D tomography model that is re-interpolated on a regular grid
 
 ## Goal
 This explains how to load a 3D seismic data set that is given in CSV format (comma separated ASCII), and plot it in paraview. The example is a shear-wave velocity model of the Alpine-Mediterranean region, described in:
diff --git a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
index fb68aeef..6ab8952e 100644
--- a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
+++ b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md
@@ -1,4 +1,4 @@
-# 3D tomography model that is given as a netCDF file
+# 2 - 3D tomography model that is given as a netCDF file
 
 ## Goal
 This explains how to load a 3D seismic data set that is given in netCDF format, and plot it in paraview. The example is a shear-wave velocity model of the Alpine-Mediterranean region, described in:
diff --git a/docs/src/man/tutorial_local_Flegrei.md b/docs/src/man/tutorial_local_Flegrei.md
index 4bbf6a02..b1303bc3 100644
--- a/docs/src/man/tutorial_local_Flegrei.md
+++ b/docs/src/man/tutorial_local_Flegrei.md
@@ -1,4 +1,4 @@
-# Km-scale volcano tutorial using cartesian coordinates
+# 13 - Campi Flegrei Volcano Tutorial using cartesian coordinates
 
 ## Goal
 
diff --git a/docs/src/man/tutorial_time_Seismicity.md b/docs/src/man/tutorial_time_Seismicity.md
index a87fe212..dadc6a2e 100644
--- a/docs/src/man/tutorial_time_Seismicity.md
+++ b/docs/src/man/tutorial_time_Seismicity.md
@@ -1,5 +1,4 @@
-# How to create a movie to show seismic distribution through time
-
+# 15 - How to create a movie that shows the temporal evolution of seismicity
 ## Goal
 
 This tutorial creates a movie of the spatial variations in seismicity using the earthquakes previously visualized at Campi Flegrei caldera. We visualize it against the travel-time model and tomography:
@@ -76,7 +75,7 @@ julia>Movie_Paraview(pvd=movie, Finalize=true)
 
 ```
 
-This tutorial has created a new *TemporalSeismicity.pvd* that can be loaded in Paraview.
+This tutorial has created a new *TemporalSeismicity.pvd* file that can be loaded in Paraview.
 
 ![Tutorial_SeismicTime_PVD](../assets/img/Tutorial_SeismicityTime_2.png)
 
diff --git a/docs/src/man/tutorials.md b/docs/src/man/tutorials.md
index 79af59d1..e99ae240 100644
--- a/docs/src/man/tutorials.md
+++ b/docs/src/man/tutorials.md
@@ -1,13 +1,4 @@
 # Tutorials
 
-The best way to learn how to use julia and GeophysicalModelGenerator to visualize your data is to look at the tutorials.
-
-1. [3D seismic tomography data on regular grid](./tutorial_load3DSeismicData.md). Demonstrates how to load 3D data that are defined on a regular longitude/latitude/depth grid.
-2. [Moho topography data](./tutorial_MohoTopo.md). Shows how to plot Moho data as 3D points in paraview, and how to fit a surface through it.
-3. [Topography](./tutorial_GMT_Topography.md). Shows how to quickly obtain the topography of any part of the world using GMT & transfer that to paraview.
-4. [Coastlines](./tutorial_Coastlines.md). Shows how to generate land/sea surfaces.
-5. [Import screenshots](./tutorial_Screenshot_To_Paraview.md). Gives examples how you can easily import screenshots from published papers and visualize them in 3D 
-6. [3D seismic tomography data on irregular grid](./tutorial_loadirregular3DSeismicData.md). Shows how to interpolate 3D seismic data, given on an irregular lon/lat grid, to a regular grid and create Paraview input from it.
-7. [Topography and geological maps](./tutorial_GMT_Topography_GeologicalMap.md). Shows how to import ETOPO1 topography, how to drape a geological map over it & transfer that to Paraview.
-8. [ISC earthquake data](./tutorial_ISC_data.md). Shows how to import earthquake data from the ISC catalogue.
-9. [Plot GPS data](./tutorial_GPS.md). Shows how to load and plot GPS data as vectors.
+The best way to learn how to use julia and GeophysicalModelGenerator to visualize your data is to look at the tutorials. We have a range of tutorials that show the functionality of `GeophysicalModelGenerator.jl` (see sidebar). They are mostly stand-alone so it doesn't matter in which order they are executed.
+The input scripts for the tutorials are available under `/tutorials`.
\ No newline at end of file
diff --git a/docs/src/man/visualise.md b/docs/src/man/visualise.md
index 7345b159..84755629 100644
--- a/docs/src/man/visualise.md
+++ b/docs/src/man/visualise.md
@@ -63,5 +63,5 @@ julia> Visualise(Data_Cart, Topography=Topo_Cart)
 ```
 
 ```@docs
-GeophysicalModelGenerator.Visualise
+Visualise
 ```
diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index c65cbc81..aa7fa64b 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -2,6 +2,7 @@ using Base: Int64, Float64, NamedTuple
 using Printf
 using Parameters        # helps setting default parameters in structures
 using SpecialFunctions: erfc
+using GeoParams
 
 # Setup_geometry
 #
@@ -10,7 +11,7 @@ using SpecialFunctions: erfc
 
 export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!,
         makeVolcTopo,
-        ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp,
+        ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
         Compute_ThermalStructure, Compute_Phase
 
@@ -36,7 +37,7 @@ Parameters
 - StrikeAngle - strike angle of slab
 - DipAngle - dip angle of slab
 - phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
-- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`
+- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`,`LithosphericTemp()`
 
 
 Examples
@@ -106,14 +107,16 @@ function AddBox!(Phase, Temp, Grid::AbstractGeneralGrid;                 # requi
                     (Yrot .>= (ylim[1] - Origin[2])) .& (Yrot .<= (ylim[2] - Origin[2])) .&
                     (Zrot .>= zbot) .& (Zrot .<= ztop)  )
 
-
     # Compute thermal structure accordingly. See routines below for different options
-    if T != nothing
-        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], T)
+    if T != nothing 
+        if isa(T,LithosphericTemp)
+            Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase)
+        end
+        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T)
     end
 
-    # Set the phase. Different routines are available for that - see below.
-    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase)
+    # Set the phase. Different routines are available for that - see below.    
+    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase)        
 
     return nothing
 end
@@ -206,7 +209,7 @@ function AddLayer!(Phase, Temp, Grid::AbstractGeneralGrid;      # required input
 
     # Compute thermal structure accordingly. See routines below for different options
     if !isnothing(T)
-        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], T)
+        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
     end
 
     # Set the phase. Different routines are available for that - see below.
@@ -273,7 +276,7 @@ function AddSphere!(Phase, Temp, Grid::AbstractGeneralGrid;      # required inpu
 
     # Compute thermal structure accordingly. See routines below for different options
     if T != nothing
-        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], T)
+        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
     end
 
     # Set the phase. Different routines are available for that - see below.
@@ -357,7 +360,7 @@ function AddEllipsoid!(Phase, Temp, Grid::AbstractGeneralGrid;      # required i
 
     # Compute thermal structure accordingly. See routines below for different options
     if T != nothing
-        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], T)
+        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T)
     end
 
     # Set the phase. Different routines are available for that - see below.
@@ -438,7 +441,7 @@ function AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid;   # required input
 
     # Compute thermal structure accordingly. See routines below for different options
     if T != nothing
-        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], T)
+        Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
     end
 
     # Set the phase. Different routines are available for that - see below.
@@ -593,7 +596,7 @@ Parameters
     T = 1000
 end
 
-function Compute_ThermalStructure(Temp, X, Y, Z, s::ConstantTemp)
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::ConstantTemp)
     Temp .= s.T
     return Temp
 end
@@ -615,7 +618,7 @@ Parameters
     Tbot = 1350
 end
 
-function Compute_ThermalStructure(Temp, X, Y, Z, s::LinearTemp)
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::LinearTemp)
     @unpack Ttop, Tbot  = s
 
     dz   = Z[end]-Z[1];
@@ -645,7 +648,7 @@ Parameters
     Adiabat = 0        # Adiabatic gradient in K/km
 end
 
-function Compute_ThermalStructure(Temp, X, Y, Z, s::HalfspaceCoolingTemp)
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp)
     @unpack Tsurface, Tmantle, Age, Adiabat  = s
 
     kappa       =   1e-6;
@@ -688,7 +691,7 @@ Parameters
     maxAge  = 60       # maximum thermal age of plate [Myrs]
 end
 
-function Compute_ThermalStructure(Temp, X, Y, Z, s::SpreadingRateTemp)
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
     @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge  = s
 
     kappa       =   1e-6;
@@ -723,7 +726,225 @@ function Compute_ThermalStructure(Temp, X, Y, Z, s::SpreadingRateTemp)
     return Temp
 end
 
+"""
+    LithosphericTemp(Tsurface=0.0, Tpot=1350.0, dTadi=0.5, 
+                        ubound="const", lbound="const, utbf = 50.0e-3, ltbf = 10.0e-3, 
+                        age = 120.0, dtfac = 0.9, nz = 201, 
+                        rheology = example_CLrheology() 
+                    )
+
+Calculates a 1D temperature profile [C] for variable thermal parameters including radiogenic heat source and 
+    linearly interpolates the temperature profile onto the box. The thermal parameters are defined in 
+    rheology and the structure of the lithosphere is define by LithosphericPhases().
+
 
+Parameters
+========
+- Tsurface  : surface temperature [C]
+- Tpot      : potential mantle temperature [C]
+- dTadi     : adiabatic gradient [K/km]
+- ubound    : Upper thermal boundary condition ["const","flux"] 
+- lbound    : Lower thermal boundary condition ["const","flux"]
+- utbf      : Upper thermal heat flux [W/m]; if ubound == "flux"
+- ltbf      : Lower thermal heat flux [W/m]; if lbound == "flux"
+- age       : age of the lithosphere [Ma]
+- dtfac     : Diffusion stability criterion to calculate T_age
+- nz        : Grid spacing for the 1D profile within the box
+- rheology  : Structure containing the thermal parameters for each phase [default example_CLrheology]
+
+"""
+@with_kw_noshow mutable struct LithosphericTemp <: AbstractThermalStructure
+    Tsurface = 0.0      # top T [C]
+    Tpot = 1350.0       # potential T [C]
+    dTadi = 0.5         # adiabatic gradient in K/km
+    ubound = "const"    # Upper thermal boundary condition
+    lbound = "const"    # lower thermal boundary condition
+    utbf = 50.0e-3      # q [W/m^2]; if ubound = "flux"
+    ltbf = 10.0e-3      # q [W/m^2]; if lbound = "flux"
+    age = 120.0         # Lithospheric age [Ma]
+    dtfac = 0.9         # Diffusion stability criterion
+    nz = 201            
+    rheology = example_CLrheology()
+end
+
+struct Thermal_parameters{A}
+    ρ::A
+    Cp::A
+    k::A
+    ρCp::A
+    H::A
+    function Thermal_parameters(ni)
+        ρ   =   zeros(ni)
+        Cp  =   zeros(ni)
+        k   =   zeros(ni)
+        ρCp =   zeros(ni)
+        H   =   zeros(ni)
+        new{typeof(ρ)}(ρ,Cp,k,ρCp,H)
+    end
+end
+
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::LithosphericTemp)
+    @unpack Tsurface, Tpot, dTadi, ubound, lbound, utbf, ltbf, age, 
+        dtfac, nz, rheology = s
+
+    # Create 1D depth profile within the box
+    z   =   LinRange(round(maximum(Z)),round(minimum(Z)),nz)    # [km]
+    z   =   @. z*1e3                                            # [m] 
+    dz  =   z[2] - z[1]                                         # Gride resolution
+
+    # Initialize 1D arrays for explicit solver
+    T       =   zeros(nz)    
+    phase   =   Int64.(zeros(nz))
+
+    # Assign phase id from Phase to 1D phase array
+    phaseid     =   (minimum(Phase):1:maximum(Phase))    
+    ztop        =   round(maximum(Z[findall(Phase .== phaseid[1])]))
+    zlayer      =   zeros(length(phaseid))
+    for i = 1:length(phaseid)
+        # Calculate layer thickness from Phase array
+        zlayer[i]   =   round(minimum(Z[findall(Phase .== phaseid[i])]))
+        zlayer[i]   =   zlayer[i]*1.0e3
+    end
+    for i = 1:length(phaseid)
+        # Assign phase ids
+        ind         =   findall((z .>= zlayer[i]) .&  (z .<= ztop))    
+        phase[ind]  .=  phaseid[i]
+        ztop        =   zlayer[i]
+    end
+
+    # Setup initial T-profile
+    Tpot        =   Tpot + 273.15                   # Potential temp [K]  
+    Tsurface    =   Tsurface + 273.15               # Surface temperature [ K ]
+    T           =   @. Tpot + abs.(z./1.0e3)*dTadi  # Initial T-profile [ K ]    
+    T[1]        =   Tsurface   
+    
+    args        = (;)
+    thermal_parameters  = Thermal_parameters(nz)   
+    
+    ## Update thermal parameters ======================================== #
+    compute_density!(thermal_parameters.ρ,rheology,phase,args)
+    compute_heatcapacity!(thermal_parameters.Cp,rheology,phase,args)
+    compute_conductivity!(thermal_parameters.k,rheology,phase,args)
+    thermal_parameters.ρCp  .=   @. thermal_parameters.Cp * thermal_parameters.ρ
+    compute_radioactive_heat!(thermal_parameters.H,rheology,phase,args)
+
+    # Thermal diffusivity [ m^2/s ]
+    κ       =   maximum(thermal_parameters.k) / 
+        minimum(thermal_parameters.ρ) / minimum(thermal_parameters.Cp)
+    ## =================================================================== #
+    ## Time stability criterion ========================================= #
+    tfac    =   60.0*60.0*24.0*365.25   # Seconds per year
+    age     =   age*1.0e6*tfac          # Age in seconds
+    dtexp   =   dz^2.0/2.0/κ            # Stability criterion for explicit
+    dt      =   dtfac*dtexp             # [s]
+    nit     =   Int64(ceil(age/dt))     # Number of iterations
+    time    =   zeros(nit)              # Time array
+    
+    for i = 1:nit
+        if i > 1
+            time[i]   =   time[i-1] + dt
+        end
+        SolveDiff1Dexplicit_vary!(
+            T,
+            thermal_parameters,
+            ubound,lbound,
+            utbf,ltbf,
+            dz,
+            dt)
+    end
+
+    interp_linear_T = linear_interpolation(-z./1.0e3, T.-273.15)      # create interpolation object
+    Temp = interp_linear_T(-Z)
+    
+    return Temp
+end
+
+function SolveDiff1Dexplicit_vary!(    
+    T,
+    thermal_parameters,
+    ubound,lbound,
+    utbf,ltbf,
+    di,
+    dt
+)    
+    nz      =   length(T)
+    T0      =   T
+
+    if ubound == "const"
+        T[1]    =   T0[1]
+    elseif ubound == "flux"
+        kB      =   (thermal_parameters.k[2] + thermal_parameters.k[1])/2.0
+        kA      =   (thermal_parameters.k[1] + thermal_parameters.k[1])/2.0
+        a       =   (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[1])
+        b       =   1 - (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[1])        
+        c       =   (dt*2.0*utbf)/(di * thermal_parameters.ρCp[1])
+        T[1]    =   a*T0[2] + b*T0[1] + c + 
+                thermal_parameters.H[1]*dt/thermal_parameters.ρCp[1]
+    end
+    if lbound == "const"
+        T[nz]   =   T0[nz]
+    elseif lbound == "flux"
+        kB      =   (thermal_parameters.k[nz] + thermal_parameters.k[nz])/2.0
+        kA      =   (thermal_parameters.k[nz] + thermal_parameters.k[nz-1])/2.0
+        a       =   (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[nz])
+        b       =   1 - (dt*(kA + kB)) / (di^2.0 * thermal_parameters.ρCp[nz])
+        c       =   -(dt*2.0*ltbf) / (di * thermal_parameters.ρCp[nz])
+        T[nz]   =   a*T0[nz-1] + b*T0[nz] + c
+    end
+
+    kAi     =   @. (thermal_parameters.k[1:end-2] + thermal_parameters.k[2:end-1])/2.0
+    kBi     =   @. (thermal_parameters.k[2:end-1] + thermal_parameters.k[3:end])/2.0
+    ai      =   @. (kBi*dt)/(di^2.0*thermal_parameters.ρCp[2:end-1])
+    bi      =   @. 1.0 - (dt*(kAi + kBi))/(di^2.0*thermal_parameters.ρCp[2:end-1])
+    ci      =   @. (kAi*dt)/(di^2.0*thermal_parameters.ρCp[2:end-1])
+    T[2:end-1]   =   @. ai*T0[3:end] + bi*T0[2:end-1] + ci*T0[1:end-2] + 
+                    thermal_parameters.H[2:end-1]*dt/thermal_parameters.ρCp[2:end-1]
+    return T    
+end
+
+function example_CLrheology(;    
+    ρM=3.0e3,           # Density [ kg/m^3 ]
+    CpM=1.0e3,          # Specific heat capacity [ J/kg/K ]
+    kM=2.3,             # Thermal conductivity [ W/m/K ]
+    HM=0.0,             # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho]
+    ρUC=2.7e3,          # Density [ kg/m^3 ]
+    CpUC=1.0e3,         # Specific heat capacity [ J/kg/K ]
+    kUC=3.0,            # Thermal conductivity [ W/m/K ]
+    HUC=617.0e-12,      # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho]
+    ρLC=2.9e3,          # Density [ kg/m^3 ]
+    CpLC=1.0e3,         # Specific heat capacity [ J/kg/K ]
+    kLC=2.0,            # Thermal conductivity [ W/m/K ]
+    HLC=43.0e-12,       # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho]
+)    
+
+    rheology = (
+        # Name              = "UpperCrust",
+        SetMaterialParams(;
+            Phase               =   1,
+            Density             =   ConstantDensity(; ρ=ρUC),
+            HeatCapacity        =   ConstantHeatCapacity(; Cp=CpUC),
+            Conductivity        =   ConstantConductivity(; k=kUC),
+            RadioactiveHeat     =   ConstantRadioactiveHeat(; H_r=HUC*ρUC),     # [H] = W/m^3
+        ),
+        # Name              = "LowerCrust",
+        SetMaterialParams(;
+            Phase               =   2,
+            Density             =   ConstantDensity(; ρ=ρLC),
+            HeatCapacity        =   ConstantHeatCapacity(; Cp=CpLC),
+            Conductivity        =   ConstantConductivity(; k=kLC),
+            RadioactiveHeat     =   ConstantRadioactiveHeat(; H_r=HLC*ρLC),     # [H] = W/m^3
+        ),
+        # Name              = "LithosphericMantle",
+        SetMaterialParams(;
+            Phase               =   3,
+            Density             =   ConstantDensity(; ρ=ρM),
+            HeatCapacity        =   ConstantHeatCapacity(; Cp=CpM),
+            Conductivity        =   ConstantConductivity(; k=kM),
+            RadioactiveHeat     =   ConstantRadioactiveHeat(; H_r=HM*ρM),       # [H] = W/m^3
+        ),
+    )
+    return rheology
+end
 
 abstract type AbstractPhaseNumber end
 
diff --git a/src/utils.jl b/src/utils.jl
index 31ad6ff9..4f8dc2cc 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -2,7 +2,7 @@
 
 export meshgrid, CrossSection, CrossSectionVolume, CrossSectionSurface, CrossSectionPoints, ExtractSubvolume, SubtractHorizontalMean
 export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap
-export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields
+export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane
 export RotateTranslateScale
 export DrapeOnTopo, LithostaticPressure!
 export FlattenCrossSection
@@ -1327,7 +1327,17 @@ function InterpolateDataOnSurface(V::GeoData, Surf::GeoData)
     return Surf_interp
 end
 
+"""    
+    Surf_interp = InterpolateTopographyOnPlane(V::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange)
 
+Interpolates a 3D data set `V` with a projection point `proj` on a plane defined by `x` and `y`, where are uniformly spaced.
+Returns the 2D array `Surf_interp`. 
+"""
+function InterpolateTopographyOnPlane(V::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange)
+    cart_grid = CartData(XYZGrid(x, y, 0))
+    tproj = ProjectCartData(cart_grid, V, proj)
+    return tproj.z.val[:, :, 1]
+end
 
 # Extracts a sub-data set using indices
 function ExtractDataSets(V::AbstractGeneralGrid, iLon, iLat, iDepth)
diff --git a/test.jld2 b/test.jld2
deleted file mode 100644
index 2386ca68..00000000
Binary files a/test.jld2 and /dev/null differ
diff --git a/test/download_GMG_temp.jld2 b/test/download_GMG_temp.jld2
deleted file mode 100644
index 8d784449..00000000
Binary files a/test/download_GMG_temp.jld2 and /dev/null differ
diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl
index 67ae3b1f..a6477d1f 100644
--- a/test/test_setup_geometry.jl
+++ b/test/test_setup_geometry.jl
@@ -1,5 +1,5 @@
 # test setting geometries in the different grid types
-using Test, GeophysicalModelGenerator
+using Test, GeophysicalModelGenerator, GeoParams
 
 
 # GeoData
@@ -96,3 +96,85 @@ ind         =   BelowSurface(Grid, Topo_cart);
 CharDim     =   GEO_units();
 Grid        =   CreateCartGrid(size=(10,20,30),x=(0.0km,10km), y=(0.0km, 10km), z=(-10.0km, 2.0km), CharDim=CharDim)
 @test Grid.Δ[2] ≈ 0.0005263157894736842
+
+
+# test 1D-explicit thermal solver for AddBox -----------
+nel         =   96
+Grid        =   CreateCartGrid(size=(nel,nel,nel),x=(-200.,200.), y=(-200.,200.), z=(-200.,0))
+Temp        =   zeros(Float64, Grid.N...);
+Phases      =   zeros(Int64,  Grid.N...);
+
+# 1) horizontally layer lithosphere; UpperCrust,LowerCrust,Mantle
+AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0), Origin=(0.0,0.0,0.0),
+    phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
+    DipAngle=0.0, T=LithosphericTemp(nz=201))
+
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 36131.638045729735
+
+# 2) inclined lithosphere; UpperCrust,LowerCrust,Mantle
+Temp    =   zeros(Float64, Grid.N...);
+Phases  =   zeros(Int64,  Grid.N...);
+
+AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0), Origin=(0.0,0.0,0.0),
+    phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
+    DipAngle=30.0, T=LithosphericTemp(nz=201))
+
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 41912.18172533137
+
+# 3) inclined lithosphere with respect to the default origin; UpperCrust,LowerCrust,Mantle
+Temp    =   zeros(Float64, Grid.N...);
+Phases  =   zeros(Int64,  Grid.N...);
+
+AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
+    phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
+    DipAngle=30.0, T=LithosphericTemp(nz=201))
+
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 41316.11499878003
+
+# 4) inclined lithosphere with only two layers
+Temp    =   zeros(Float64, Grid.N...);
+Phases  =   zeros(Int64,  Grid.N...);
+
+ρM=3.0e3            # Density [ kg/m^3 ]
+CpM=1.0e3           # Specific heat capacity [ J/kg/K ]
+kM=2.3              # Thermal conductivity [ W/m/K ]
+HM=0.0              # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho]
+ρUC=2.7e3           # Density [ kg/m^3 ]
+CpUC=1.0e3          # Specific heat capacity [ J/kg/K ]
+kUC=3.0             # Thermal conductivity [ W/m/K ]
+HUC=617.0e-12       # Radiogenic heat source per mass [H] = W/kg; [H] = [Q/rho]
+
+rheology = (
+        # Name              = "UpperCrust",
+        SetMaterialParams(;
+            Phase               =   1,
+            Density             =   ConstantDensity(; ρ=ρUC),
+            HeatCapacity        =   ConstantHeatCapacity(; Cp=CpUC),
+            Conductivity        =   ConstantConductivity(; k=kUC),
+            RadioactiveHeat     =   ConstantRadioactiveHeat(; H_r=HUC*ρUC),     # [H] = W/m^3
+        ),
+        # Name              = "LithosphericMantle",
+        SetMaterialParams(;
+            Phase               =   2,
+            Density             =   ConstantDensity(; ρ=ρM),
+            HeatCapacity        =   ConstantHeatCapacity(; Cp=CpM),
+            Conductivity        =   ConstantConductivity(; k=kM),
+            RadioactiveHeat     =   ConstantRadioactiveHeat(; H_r=HM*ρM),       # [H] = W/m^3
+        ),
+    );
+
+AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
+    phase=LithosphericPhases(Layers=[20 80], Phases = [1 2], Tlab=nothing), 
+    DipAngle=30.0, T=LithosphericTemp(rheology=rheology,nz=201))
+
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 40513.969831615716
+
+# 5) using flux lower boundary conditions
+Temp    =   zeros(Float64, Grid.N...);
+Phases  =   zeros(Int64,  Grid.N...);
+
+AddBox!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
+    phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing), 
+    DipAngle=30.0, T=LithosphericTemp(lbound="flux",nz=201))
+
+@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 37359.648604722104
\ No newline at end of file
diff --git a/tutorial/Alps_VpModel_Zhao_etal_JGR2016.jl b/tutorials/Alps_VpModel_Zhao_etal_JGR2016.jl
similarity index 100%
rename from tutorial/Alps_VpModel_Zhao_etal_JGR2016.jl
rename to tutorials/Alps_VpModel_Zhao_etal_JGR2016.jl
diff --git a/docs/src/scripts/LaPalma_example.jl b/tutorials/LaPalma_example.jl
similarity index 93%
rename from docs/src/scripts/LaPalma_example.jl
rename to tutorials/LaPalma_example.jl
index 4ce8b412..81e50af4 100644
--- a/docs/src/scripts/LaPalma_example.jl
+++ b/tutorials/LaPalma_example.jl
@@ -1,7 +1,8 @@
-# # Create a 3D LaMEM setup for La Palma
+# # Create a Cartesian Model Setup for La Palma
 #
 # ## Aim
-# In this tutorial, your will learn how to use real data to create a geodynamic model setup with LaMEM. We will use the data of La Palma, which is a volcanic island that started erupting in mid september 2021.
+# In this tutorial, your will learn how to use real data to create a geodynamic model setup with LaMEM. 
+# We will use the data of the Cumbre Viejo eruption in La Palma, which is a volcanic island that erupted from mid september 2021 - december 2021.
 # LaMEM is a cartesian geodynamic model, which implies that you will have to transfer the data from `GeoData` to `CartData`.
 #
 #
@@ -10,9 +11,8 @@
 #  1) Topography
 #  2) Earthquake locations
 
-# We start with loading the required packages
-using GeophysicalModelGenerator, JLD2
-using GMT
+# We start with loading the required packages, which includes `GMT` to download topography (an optional dependency for `GeophysicalModelGenerator`)
+using GeophysicalModelGenerator, GMT
 
 # In case you managed to install GMT on your machine, you can automatically download the topography with
 Topo = ImportTopo(lon = [-18.7, -17.1], lat=[28.0, 29.2], file="@earth_relief_03s.grd")
@@ -57,6 +57,9 @@ Topo_LaMEM = ProjectCartData(Topo_LaMEM, Topo, proj)
 Write_Paraview(EQ_cart,"EQ_cart",PointsData=true)
 Write_Paraview(Topo_LaMEM,"Topo_LaMEM")
 
+
+#= 
+
 # ## 3. Create LaMEM setup
 #
 # In a next step, we need to read the LaMEM input file of the simulation. In particular, this will read the lateral dimensions of the grid and
@@ -123,3 +126,5 @@ Save_LaMEMMarkersParallel(Model3D)
 
 #src Note: The markdown page is generated using:
 #src Literate.markdown("LaPalma_example.jl",outputdir="../man/", execute=true)
+
+=#
\ No newline at end of file
diff --git a/tutorial/Lippitsch_Screenshots.jl b/tutorials/Lippitsch_Screenshots.jl
similarity index 100%
rename from tutorial/Lippitsch_Screenshots.jl
rename to tutorials/Lippitsch_Screenshots.jl
diff --git a/tutorial/MeRe_ElSharkawy.jl b/tutorials/MeRe_ElSharkawy.jl
similarity index 100%
rename from tutorial/MeRe_ElSharkawy.jl
rename to tutorials/MeRe_ElSharkawy.jl
diff --git a/tutorial/MohoTopo_Spada.jl b/tutorials/MohoTopo_Spada.jl
similarity index 100%
rename from tutorial/MohoTopo_Spada.jl
rename to tutorials/MohoTopo_Spada.jl
diff --git a/tutorial/Tutorial_AlpineData.jl b/tutorials/Tutorial_AlpineData.jl
similarity index 100%
rename from tutorial/Tutorial_AlpineData.jl
rename to tutorials/Tutorial_AlpineData.jl
diff --git a/tutorial/Tutorial_CSV.jl b/tutorials/Tutorial_CSV.jl
similarity index 100%
rename from tutorial/Tutorial_CSV.jl
rename to tutorials/Tutorial_CSV.jl
diff --git a/tutorial/Tutorial_Flegrei.jl b/tutorials/Tutorial_Flegrei.jl
similarity index 100%
rename from tutorial/Tutorial_Flegrei.jl
rename to tutorials/Tutorial_Flegrei.jl
diff --git a/tutorial/Tutorial_GPS.jl b/tutorials/Tutorial_GPS.jl
similarity index 100%
rename from tutorial/Tutorial_GPS.jl
rename to tutorials/Tutorial_GPS.jl
diff --git a/tutorial/Tutorial_SeismicityTime.jl b/tutorials/Tutorial_SeismicityTime.jl
similarity index 100%
rename from tutorial/Tutorial_SeismicityTime.jl
rename to tutorials/Tutorial_SeismicityTime.jl
diff --git a/tutorial/Tutorial_TopoGeological_PNG.jl b/tutorials/Tutorial_TopoGeological_PNG.jl
similarity index 100%
rename from tutorial/Tutorial_TopoGeological_PNG.jl
rename to tutorials/Tutorial_TopoGeological_PNG.jl
diff --git a/tutorial/Tutorial_TopoGeological_SHP.jl b/tutorials/Tutorial_TopoGeological_SHP.jl
similarity index 100%
rename from tutorial/Tutorial_TopoGeological_SHP.jl
rename to tutorials/Tutorial_TopoGeological_SHP.jl
diff --git a/docs/src/scripts/Tutorial_Votemaps.jl b/tutorials/Tutorial_Votemaps.jl
similarity index 99%
rename from docs/src/scripts/Tutorial_Votemaps.jl
rename to tutorials/Tutorial_Votemaps.jl
index 597f26cc..1aebf5d6 100644
--- a/docs/src/scripts/Tutorial_Votemaps.jl
+++ b/tutorials/Tutorial_Votemaps.jl
@@ -1,4 +1,4 @@
-# # Votemaps
+# # 12 - Votemaps
 #
 # ## Aim
 # In this tutorial, your will learn how to create Votemaps that compare different tomographic models and look for similarities between different models.
diff --git a/tutorial/Tutorial_netCDF.jl b/tutorials/Tutorial_netCDF.jl
similarity index 100%
rename from tutorial/Tutorial_netCDF.jl
rename to tutorials/Tutorial_netCDF.jl