diff --git a/Khefalonia2D.jl b/Khefalonia2D.jl
index 72e09402..5b734be9 100644
--- a/Khefalonia2D.jl
+++ b/Khefalonia2D.jl
@@ -131,7 +131,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
 
     # PT coefficients for thermal diffusion
     pt_thermal       = PTThermalCoeffs(
-        backend, rheology, phase_ratios, args0, dt, ni, di, li; ϵ=1e-5, CFL=0.5 / √2
+        backend, rheology, phase_ratios, args0, dt, ni, di, li; ϵ=1e-5, CFL=1e-2 / √2
     )
 
     # Boundary conditions
@@ -224,7 +224,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
                     iterMax          = 50e3,
                     nout             = 1e3,
                     viscosity_cutoff = viscosity_cutoff,
-                    free_surface     = false,
+                    free_surface     = true,
                     viscosity_relaxation = 1e-2
                 )
             );
@@ -285,6 +285,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
                 data_v = (;
                     T   = Array(T_buffer),
                     τII = Array(stokes.τ.II),
+                    τxx = Array(stokes.τ.xx),
+                    τyy = Array(stokes.τ.yy),
                     εII = Array(stokes.ε.II),
                     Vx  = Array(Vx_v),
                     Vy  = Array(Vy_v),
diff --git a/KhefaloniaRheology.jl b/KhefaloniaRheology.jl
index 44956b53..9e2dfc49 100644
--- a/KhefaloniaRheology.jl
+++ b/KhefaloniaRheology.jl
@@ -21,7 +21,7 @@ function init_rheologies(; is_plastic = true)
     # β_sublithospheric_mantle    = inv(get_Kb(el_sublithospheric_mantle))
 
     # Physical properties using GeoParams ----------------
-    η_reg     = 1e16
+    η_reg     = 1e17
     cohesion  = 15e6
     friction  = 20
     pl_crust  = DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
diff --git a/KhefaloniaSetup.jl b/KhefaloniaSetup.jl
index ba67d529..0e76574b 100644
--- a/KhefaloniaSetup.jl
+++ b/KhefaloniaSetup.jl
@@ -83,7 +83,7 @@ function Khefalonia__setup(N)
     # @. Temp               = max(Temp, 20)
     Grid = addfield(Grid,(; Phases, Temp))
 
-    li     = abs(last(x)-first(x)),  abs(last(z)-first(z))
+    li     = (abs(last(x)-first(x)),  abs(last(z)-first(z))) .* 1e3
     origin = x[1], z[1]
 
     ph      = Phases