diff --git a/TurbulenceModels/wallFunctions b/TurbulenceModels/wallFunctions
index b12e4d362..a1e33e27a 160000
--- a/TurbulenceModels/wallFunctions
+++ b/TurbulenceModels/wallFunctions
@@ -1 +1 @@
-Subproject commit b12e4d3621915cf12b2636a8f18ceebd82417e7b
+Subproject commit a1e33e27ae86196294c692f5f207bfe33023716e
diff --git a/docker/docker-stacks b/docker/docker-stacks
index 9abc77b95..c2159f411 160000
--- a/docker/docker-stacks
+++ b/docker/docker-stacks
@@ -1 +1 @@
-Subproject commit 9abc77b951c45a293a7d9b747163ff682e38bed5
+Subproject commit c2159f4118f39e7201282628131dad3b25d40ba6
diff --git a/tutorials/Py/plot_depth2DPipelineScour.py b/tutorials/Py/plot_depth2DPipelineScour.py
index e201fc389..659e26d24 100644
--- a/tutorials/Py/plot_depth2DPipelineScour.py
+++ b/tutorials/Py/plot_depth2DPipelineScour.py
@@ -2,6 +2,7 @@
 from pylab import mlab
 import matplotlib.pyplot as plt
 import fluidfoam
+from scipy.interpolate import griddata
 
 
 # Function reading the mesh dimensions
@@ -12,7 +13,7 @@ def readxy(case):
 
 # Function returning the bed profile
 def max_depth(case, t0, tfinal, dt, xi, yi):
-    ybed = np.zeros(len(xi))
+    ybed = np.zeros(ngridx)
     arr_size = int((tfinal - t0) / dt + 0.01)
     yscour = np.zeros(arr_size)
     time = np.zeros(arr_size)
@@ -26,10 +27,10 @@ def max_depth(case, t0, tfinal, dt, xi, yi):
             timename = str(int(t)) + "/"
         else:
             timename = str(t) + "/"
-        a = fluidfoam.readscalar(case, timename, "alpha_a")
-        ai = mlab.griddata(x, y, a, xi, yi, interp="linear")
+        a = fluidfoam.readscalar(case, timename, "alpha.a")
+        ai = griddata((x, y), a, (xi, yi), method="linear")
         for j in range(ngridx):
-            ybed[j] = yi[np.max((np.where(ai[:, j] > 0.5)))]
+            ybed[j] = np.max(yi[(np.where(ai[:, j] > 0.5))])
         yscour[i] = np.min(ybed[:])
         time[i] = t
     return yscour, time
@@ -51,13 +52,14 @@ def max_depth(case, t0, tfinal, dt, xi, yi):
 # Interpolation grid
 xi = np.linspace(xinterpmin, xinterpmax, ngridx)
 yi = np.linspace(yinterpmin, yinterpmax, ngridy)
+xinterp, yinterp = np.meshgrid(xi, yi)
 
 # Maximum depth calculation
 x, y = readxy(case)
 t0 = 0
 tfinal = 30
 dt = 0.5
-max_depth, time = max_depth(case, t0, tfinal, dt, xi, yi)
+max_depth, time = max_depth(case, t0, tfinal, dt, xinterp, yinterp)
 
 # Experimental data collection
 expe = "DATA/Mao_depth_expe.txt"
diff --git a/tutorials/Py/plot_profiles2DPipelineScour.py b/tutorials/Py/plot_profiles2DPipelineScour.py
index 18dcdda17..a9ea85fc5 100644
--- a/tutorials/Py/plot_profiles2DPipelineScour.py
+++ b/tutorials/Py/plot_profiles2DPipelineScour.py
@@ -19,7 +19,7 @@ def depth(sol, t, x, y, xi, yi):
         timename = str(int(t)) + "/"
     else:
         timename = str(t) + "/"
-    alpha = fluidfoam.readscalar(sol, timename, "alpha_a")
+    alpha = fluidfoam.readscalar(sol, timename, "alpha.a")
     alphai = griddata((x, y), alpha, (xi, yi))
     for j in range(Nx - 1):
         tab = np.where(alphai[:, j + 1] > 0.5)
diff --git a/tutorials/RAS/2DPipelineScour/constant/kineticTheoryProperties b/tutorials/RAS/2DPipelineScour/constant/kineticTheoryProperties
index d1153bada..1385a1fd1 100644
--- a/tutorials/RAS/2DPipelineScour/constant/kineticTheoryProperties
+++ b/tutorials/RAS/2DPipelineScour/constant/kineticTheoryProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-kineticTheory   on;
+kineticTheory   off;
 
 limitProduction off;