diff --git a/notebooks/Rivulet_analysis.jl b/notebooks/Rivulet_analysis.jl
index f5d97bc..742c33f 100644
--- a/notebooks/Rivulet_analysis.jl
+++ b/notebooks/Rivulet_analysis.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.19.32
+# v0.19.36
 
 using Markdown
 using InteractiveUtils
@@ -515,7 +515,7 @@ function t0_data()
 	# Loop through initial conditions
 	for angle in [2/9, 1/6, 1/9, 1/18]
 		for R in [150, 160, 180, 200]
-			for rr in [20, 30, 40, 80]
+			for rr in [20, 30, 40, 60, 80, 100]
 				# Create initial condition
 				h = torus(512, 512, rr, R, angle, (256, 256))
 				# measure relevant parameter
@@ -564,7 +564,7 @@ end
 
 Generates video files in gif formate for the supplied set of `data`.
 """
-function data2gif(data; prefix="")
+function data2gif(data; prefix="", arrested=false)
 	kbtDict = Dict(0.0 => "kbt_off", 1.0e-6 => "kbt_on")
 	maxtimestep = 2500000
 	for i in eachindex(data)
@@ -572,7 +572,7 @@ function data2gif(data; prefix="")
 		if isfile(filename)
 			println("There is alread a file called $(prefix)ang_$(data[i][8])_R_$(data[i][1])_rr_$(data[i][2])_$(kbtDict[data[i][3]]).gif  in the assets folder")
 		else
-			h = read_data(R=data[i][1], r=data[i][2], kbT=data[i][3], month=data[i][4], day=data[i][5], hour=data[i][6], minute=data[i][7], θ=data[i][8], nm=32)
+			h = read_data(R=data[i][1], r=data[i][2], kbT=data[i][3], month=data[i][4], day=data[i][5], hour=data[i][6], minute=data[i][7], θ=data[i][8], nm=32, arrested=arrested)
 			println("R=$(data[i][1]) with rr=$(data[i][2]) and kbt=$(data[i][3])")
 			if data[i][3] == 0.0
 				do_gif(h, "$(prefix)ang_$(data[i][8])_R_$(data[i][1])_rr_$(data[i][2])_$(kbtDict[data[i][3]])", timeMax=maxtimestep)
@@ -731,12 +731,20 @@ data_arrested = [
 	(160, 20, 0.0,    11,25, 14, 19, 30), 	#
 	(160, 30, 0.0,    11,25, 15, 42, 30), 	#
 	(160, 40, 0.0,    11,25, 17,  4, 30), 	#
+	(160, 60, 0.0,    12,23,  0, 32, 30), 	#
+	(160, 80, 0.0,    12,23,  1, 53, 30), 	#
 	(180, 20, 0.0,    11,25, 18, 27, 30), 	#
 	(180, 30, 0.0,    11,25, 19, 48, 30), 	#
 	(180, 40, 0.0,    11,25, 21,  7, 30), 	#
+	(180, 60, 0.0,    12,23,  4, 34, 30), 	#
+	(180, 80, 0.0,    12,23,  5, 55, 30), 	#
+	(180, 100,0.0,    12,23,  7, 15, 30), 	#
 	(200, 20, 0.0,    11,25, 22, 23, 30), 	#
 	(200, 30, 0.0,    11,25, 23, 40, 30), 	#
 	(200, 40, 0.0,    11,26,  1,  2, 30), 	#
+	(200, 60, 0.0,    12,23,  8, 36, 30), 	#
+	(200, 80, 0.0,    12,23,  9, 57, 30), 	#
+	(200, 100,0.0,    12,23, 11, 19, 30), 	#
 	(160, 20, 1.0e-6, 11,26,  2, 24, 30), 	#
 	(160, 30, 1.0e-6, 11,26,  3, 45, 30), 	#
 	(160, 40, 1.0e-6, 11,26,  5,  7, 30), 	#
@@ -749,12 +757,21 @@ data_arrested = [
 	(160, 20, 0.0,    11,27, 14, 51, 20), 	#
 	(160, 30, 0.0,    11,27, 16, 12, 20), 	#
 	(160, 40, 0.0,    11,27, 17, 33, 20), 	#
+	(160, 60, 0.0,    12,22, 12, 25, 20), 	#
+	(160, 80, 0.0,    12,22, 13, 44, 20), 	#
+	(160, 100,0.0,    12,22, 15,  5, 20), 	#
 	(180, 20, 0.0,    11,27, 18, 54, 20), 	#
 	(180, 30, 0.0,    11,27, 20, 15, 20), 	#
 	(180, 40, 0.0,    11,27, 21, 36, 20), 	#
+	(180, 60, 0.0,    12,22, 16, 25, 20), 	#
+	(180, 80, 0.0,    12,22, 17, 47, 20), 	#
+	(180,100, 0.0,    12,22, 19,  8, 20), 	#
 	(200, 20, 0.0,    11,27, 22, 57, 20), 	#
 	(200, 30, 0.0,    11,28,  0, 18, 20), 	#
 	(200, 40, 0.0,    11,28,  1, 39, 20), 	#
+	(200, 60, 0.0,    12,22, 20, 29, 20), 	#
+	(200, 80, 0.0,    12,22, 21, 50, 20), 	#
+	(200,100, 0.0,    12,22, 23, 11, 20), 	#
 	(160, 20, 1.0e-6, 11,28,  3,  0, 20), 	#
 	(160, 30, 1.0e-6, 11,28,  4, 22, 20), 	#
 	(160, 40, 1.0e-6, 11,28,  5, 42, 20), 	#
@@ -767,12 +784,21 @@ data_arrested = [
 	(160, 20, 0.0,    11,26, 14, 38, 10), 	#
 	(160, 30, 0.0,    11,26, 15, 59, 10), 	#
 	(160, 40, 0.0,    11,26, 17, 21, 10), 	#
+	(160, 60, 0.0,    12,23, 12, 40, 10), 	#
+	(160, 80, 0.0,    12,23, 14,  1, 10), 	#
+	(160, 100,0.0,    12,23, 15, 22, 10), 	#
 	(180, 20, 0.0,    11,26, 18, 42, 10), 	#
 	(180, 30, 0.0,    11,26, 20,  3, 10), 	#
 	(180, 40, 0.0,    11,26, 21, 24, 10), 	#
+	(180, 60, 0.0,    12,23, 16, 42, 10), 	#
+	(180, 80, 0.0,    12,23, 18,  2, 10), 	#
+	(180, 100,0.0,    12,23, 19, 21, 10), 	#
 	(200, 20, 0.0,    11,26, 22, 45, 10), 	#
 	(200, 30, 0.0,    11,27,  0,  5, 10), 	#
 	(200, 40, 0.0,    11,27,  1, 26, 10), 	#
+	(200, 60, 0.0,    12,23, 20, 35, 10), 	#
+	(200, 80, 0.0,    12,23, 21, 51, 10), 	#
+	(200, 100,0.0,    12,23, 23, 10, 10), 	#
 	(160, 20, 1.0e-6, 11,27,  2, 47, 10), 	#
 	(160, 30, 1.0e-6, 11,27,  4,  8, 10), 	#
 	(160, 40, 1.0e-6, 11,27,  5, 25, 10), 	#
@@ -1159,6 +1185,9 @@ for sims in [(data, ""), (data_arrested, "arr_")]
 	data2gif(sims[1], prefix=sims[2])
 end
 
+# ╔═╡ 0b5d27c1-aa9c-4fe6-a592-52eff9f5ee20
+data2gif(data_arrested, prefix="arr_", arrested=true)
+
 # ╔═╡ ecba3acb-6bc1-4722-9cee-a388a2442fae
 # ╠═╡ disabled = true
 #=╠═╡
@@ -3203,14 +3232,14 @@ version = "1.4.1+1"
 # ╠═f268582b-0756-41cf-910d-7a57b698451d
 # ╟─1b26468c-b4f7-4252-b891-4f95bc04c869
 # ╟─6d3c1725-75fa-412e-9b30-8f8df4e7874b
-# ╟─0acf9712-b27c-40c8-9bec-64d6389ce2c4
+# ╠═0acf9712-b27c-40c8-9bec-64d6389ce2c4
 # ╟─eadae383-6b5b-4e4e-80b9-5eb2fc4a5ead
 # ╟─789e9f0e-863a-4cd5-8f99-f830120e8960
 # ╟─2df8c833-7ca7-4d7a-ade5-0df083a013a1
 # ╟─81d255ea-1ab3-4635-ab4c-66100a820b28
 # ╟─6e82547e-c935-4a3e-b736-a0dae07bfb50
 # ╟─9da027de-9ee2-487c-b978-cbfd77e35fef
-# ╟─f25c4971-572f-41e2-be87-ad513c86e737
+# ╠═f25c4971-572f-41e2-be87-ad513c86e737
 # ╟─04e344a3-3d5b-449e-9222-481df24015c7
 # ╠═974c334e-38fb-436e-842b-bb016854d136
 # ╟─0dd77f66-ec27-4cdf-81e8-bcecfbcfcf29
@@ -3231,13 +3260,13 @@ version = "1.4.1+1"
 # ╟─70da13b0-6111-4d5d-a6f5-49fcc0499738
 # ╟─5baa1023-db81-4374-913d-1e88bacdb2ac
 # ╟─b3ae647b-1de9-4f56-b786-8719705c1e09
-# ╟─37756334-7859-499d-b355-658349aa1805
+# ╠═37756334-7859-499d-b355-658349aa1805
 # ╟─063757cb-b822-44e3-8a2b-57808c6f30cf
 # ╟─4fb1d7ad-47f2-4adf-a2ba-0ecc0fc8eeb0
 # ╟─41aee571-9016-4759-859a-c99eb143a410
 # ╟─13ce2bea-889f-4727-a126-71a5006a86ab
 # ╟─c9572357-8d97-47a7-914a-91c0b452eb6b
-# ╟─846ebcbe-34d6-48a4-bc23-cbd04bacf526
+# ╠═846ebcbe-34d6-48a4-bc23-cbd04bacf526
 # ╟─0f204a06-71b2-438a-bb49-4af8ebda0001
 # ╟─d5152b67-bc1d-4cc0-b73e-90d79dbadcb4
 # ╟─ab5b4c7c-ae24-4aae-a528-1dc427a7f1f1
@@ -3268,7 +3297,8 @@ version = "1.4.1+1"
 # ╠═512e1060-eee5-4374-966c-02d7fb62f303
 # ╠═61474944-c347-448a-beb9-aa2e4ef6331e
 # ╠═df08506c-f66d-430a-b235-4c9dfb80d414
-# ╟─44198158-5647-4074-95a6-96627894fadd
+# ╠═44198158-5647-4074-95a6-96627894fadd
+# ╠═0b5d27c1-aa9c-4fe6-a592-52eff9f5ee20
 # ╠═ecba3acb-6bc1-4722-9cee-a388a2442fae
 # ╠═87c627a8-2c54-44ff-aa66-d9b5c379f646
 # ╟─c848d2cf-5d36-4437-b53a-e278150e75ef
diff --git a/scripts/rivulet_stability.jl b/scripts/Rivulet_stability.jl
similarity index 68%
rename from scripts/rivulet_stability.jl
rename to scripts/Rivulet_stability.jl
index fe2d574..3c7a4c7 100644
--- a/scripts/rivulet_stability.jl
+++ b/scripts/Rivulet_stability.jl
@@ -29,6 +29,7 @@ function rivulet_run(
     device::String;
     shape = :ring, 
     arrested = false,
+    gradient = false,
     R = 150,
     rr = 100,
     ϵ = 0.01,
@@ -55,13 +56,24 @@ function rivulet_run(
                 theta[i] = 1/3
             end
         end
-    end
+    elseif gradient
+        # Wettability gradient that radial decreases the contact angle
+        theta = zeros(sys.Lx, sys.Ly)
+        dist = zeros(sys.Lx, sys.Ly)
+        for i in 1:sys.Lx
+            for j in 1:sys.Ly
+                dist[i,j] = round(Int, sqrt((i - Lx÷2)^2 + (j - Ly÷2)^2))
+            end
+        end
+        theta .= π/(6R) .* dist .+ π/18
+        theta[dist .> R] .= 2π/9 
+    end 
         # Push it to the desired device
     if device == "CPU"
         state.height .= h
     elseif device == "GPU"
         CUDA.copyto!(state.height, h)
-        if arrested
+        if arrested || gradient
             pinned = CUDA.zeros(Float64, sys.Lx, sys.Ly)
             CUDA.copyto!(pinned, theta)
         end
@@ -79,7 +91,7 @@ function rivulet_run(
             end
         end
         
-        if arrested
+        if arrested || gradient
             Swalbe.filmpressure!(state.pressure, state.height, state.dgrad, sys.param.γ, pinned, sys.param.n, sys.param.m, sys.param.hmin, sys.param.hcrit)
         else
             Swalbe.filmpressure!(state, sys)
@@ -108,13 +120,16 @@ end
 timeInterval = 25000
 
 # Make a parameter sweep
-for ang in [1/9, 2/9, 1/6, 1/18] # 
-    for gamma in [(0.005, "05"), (0.02, "20")] # , 1e-6
-        sys = Swalbe.SysConst(512, 512, Swalbe.Taumucs(Tmax=2500000, kbt=0.0, γ=gamma[1], n=3, m=2, θ=ang))
-        for outerRad in [160, 180, 200]
-            for innerRad in [20, 30, 40]
-                # Run the simulation
-                fluid = rivulet_run(sys, "GPU", R=outerRad, rr=innerRad, arrested=false, dump=timeInterval)
+for ang in [2/9] # 1/9, 1/6,  
+    for deltas in [1.0] # 0.5, 2.5
+        sys = Swalbe.SysConst(512, 512, Swalbe.Taumucs(Tmax=2500000, δ=deltas , n=3, m=2, θ=ang))
+        for outerRad in [180]# [160, 180, 200]
+            for innerRad in [20]# [60, 80, 100]
+            # Run the simulation
+                arr = false #true
+                grad = true #true
+                slips = false
+                fluid = rivulet_run(sys, "GPU", R=outerRad, rr=innerRad, arrested=arr, dump=timeInterval, gradient=grad)
                 df_fluid = Dict()
                 nSnapshots = sys.param.Tmax ÷ timeInterval
                 for t in 1:nSnapshots
@@ -123,12 +138,20 @@ for ang in [1/9, 2/9, 1/6, 1/18] #
                 end
                 println("Saving rivulet snapshots for R=$(outerRad) and r=$(innerRad) to disk")
                 save_ang = Int(round(rad2deg(π*sys.param.θ)))
-                file_name = "data/Rivulets/gamma$(gamma[2])_height_R_$(outerRad)_r_$(innerRad)_ang_$(save_ang)_kbt_$(sys.param.kbt)_nm_$(sys.param.n)-$(sys.param.m)_runDate_$(year(today()))$(month(today()))$(day(today()))$(hour(now()))$(minute(now())).jld2"
+                if arr
+                    file_name = "data/Rivulets/arrested_height_R_$(outerRad)_r_$(innerRad)_ang_$(save_ang)_kbt_$(sys.param.kbt)_nm_$(sys.param.n)-$(sys.param.m)_runDate_$(year(today()))$(month(today()))$(day(today()))$(hour(now()))$(minute(now())).jld2"
+                elseif slips
+                    file_name = "data/Rivulets/slip_$(Int(10*deltas))_height_R_$(outerRad)_r_$(innerRad)_ang_$(save_ang)_kbt_$(sys.param.kbt)_nm_$(sys.param.n)-$(sys.param.m)_runDate_$(year(today()))$(month(today()))$(day(today()))$(hour(now()))$(minute(now())).jld2"
+                elseif grad
+                    file_name = "data/Rivulets/wet_grad_lin_1040_height_R_$(outerRad)_r_$(innerRad)_ang_$(save_ang)_kbt_$(sys.param.kbt)_nm_$(sys.param.n)-$(sys.param.m)_runDate_$(year(today()))$(month(today()))$(day(today()))$(hour(now()))$(minute(now())).jld2"
+                else
+                    file_name = "data/Rivulets/height_R_$(outerRad)_r_$(innerRad)_ang_$(save_ang)_kbt_$(sys.param.kbt)_nm_$(sys.param.n)-$(sys.param.m)_runDate_$(year(today()))$(month(today()))$(day(today()))$(hour(now()))$(minute(now())).jld2"
+                end
                 save(file_name, df_fluid)
                 CUDA.reclaim()
                 fluid .= 0.0
                 df_fluid = Dict()
-                println("Done with $(ang) $(gamma[1]) $(outerRad) $(innerRad)")
+                println("Done with $(ang) $(outerRad) $(innerRad)")
             end
         end
     end