diff --git a/.gitignore b/.gitignore
index 5a795cc..cec82ae 100644
--- a/.gitignore
+++ b/.gitignore
@@ -42,4 +42,4 @@ __pycache__
 .pytest_cache
 .ipynb_checkpoints
 .vscode
-.gitconfig.local
+.gitconfig.local
\ No newline at end of file
diff --git a/code/Project.toml b/code/Project.toml
old mode 100644
new mode 100755
index 9bc3df8..779c848
--- a/code/Project.toml
+++ b/code/Project.toml
@@ -2,9 +2,14 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
+FreezeCurves = "71e4ad71-e4f2-45a3-aa0a-91ffaa9676be"
+IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
diff --git a/code/sensitivity/1DHeatEquation_with_Phase_Change.jl b/code/sensitivity/1DHeatEquation_with_Phase_Change.jl
new file mode 100755
index 0000000..aee2631
--- /dev/null
+++ b/code/sensitivity/1DHeatEquation_with_Phase_Change.jl
@@ -0,0 +1,126 @@
+using ComponentArrays
+using FreezeCurves
+using IfElse
+using LinearAlgebra
+using OrdinaryDiffEq
+using SciMLSensitivity
+using Statistics
+using Zygote
+
+# plotting
+import Plots
+
+const L = 3.34e8 # volumetric latent heat of fusion of water [J/m^3]
+
+const default_p = ComponentVector(
+    θwi=0.5,  # total water+ice content; 1.0 would correspond to a block of pure ice/water
+    k_s=2.5,  # thermal conductivity of solid material [W/(m.K)]
+    k_w=0.57, # thermal conductivity of unfrozen water [W/(m.K)]
+    k_i=2.2,  # thermal conductivity of ice [W/(m.K)]
+    c_s=2.5e6, # heat capacity of solid material [J/(m^3.K)]
+    c_w=4.2e6, # heat capacity of unfrozen water [J/(m^3.K)]
+    c_i=1.9e6, # heat capacity of ice [J/(m^3.K)]
+    jH_lb=0.05, # lower boundary flux [W/m^2]
+)
+
+struct Grid1D{TA}
+    edges::TA
+    cells::TA
+    thick::TA
+    dists::TA
+end
+function Grid1D(grid_edges::AbstractVector)
+    cells = (grid_edges[1:end-1] .+ grid_edges[2:end]) ./ 2
+    thick = diff(grid_edges)
+    dists = diff(cells)
+    return Grid1D(grid_edges, cells, thick, dists)
+end
+
+function steady_state_init(grid::Grid1D, p; Tsurf=-10.0)
+    T0 = Tsurf .+ p.jH_lb / p.k_s.*grid.cells
+    H0 = enthalpy.(T0, Ref(p))
+    return T0, H0
+end
+
+function heatcapacity(θw, p)
+    c = (1 - p.θwi)*p.c_s + θw*p.c_w + (1-θw)*p.c_i
+    return c
+end
+
+function enthalpy(T, p)
+    if T < zero(T)
+        T*heatcapacity(zero(p.θwi), p)
+    else
+        T*heatcapacity(p.θwi, p) + L*p.θwi
+    end
+end
+
+function enthalpyinv(H, p)
+    θwi = p.θwi
+    θw, Lθ = FreezeCurves.freewater(H, θwi, L)
+    C = heatcapacity(θw, p)
+    T_f = H / C
+    T_t = (H - Lθ) / C
+    T = IfElse.ifelse(
+        H < zero(θwi),
+        # Case 1: H < 0 -> frozen
+        T_f,
+        # Case 2: H >= 0
+        IfElse.ifelse(
+            H >= Lθ,
+            # Case 2a: H >= Lθ -> thawed
+            T_t,
+            # Case 2b: 0 <= H < Lθ -> phase change
+            zero(T_f)
+        )
+    )
+end
+
+function heateq_with_phase_change(
+    grid::Grid1D,
+    T_ub::F=(u,p,t) -> zero(eltype(u))
+) where {F}
+    function f(u, p, t)
+        Hinv(H) = enthalpyinv(H, p)
+        Hinv_res = Hinv.(u)
+        T = map(first, Hinv_res)
+        θw = map(last, Hinv_res)
+        θwi = p.θwi
+        k = (1-θwi)*p.k_s .+ θw*p.k_w .+ (1.0.-θw)*p.k_i
+        jH = k[1:end-1].*(T[1:end-1] .- T[2:end]) ./ grid.dists
+        jH_lb = p.jH_lb
+        jH_ub = -2*(T[1] - T_ub(u,p,t)) / grid.thick[1]
+        dH = vcat(
+            (jH_ub - jH[1]) / grid.thick[1],
+            (jH[1:end-1] .- jH[2:end]) ./ grid.thick[2:end-1],
+            (jH[end] - jH_lb) / grid.thick[end]
+        )
+        return zero(u) .+ dH
+    end
+    return f
+end
+
+# upper boundary temperature
+T_ub(u,p,t) = 10*sin(2π*t/(24*3600))
+
+p = default_p
+grid_edges = vcat(0.0:0.05:5.0, 5.1:0.1:20.0, 20.5:0.5:50.0, 51.0:1.0:100.0)
+grid = Grid1D(grid_edges)
+T0, H0 = steady_state_init(grid, p; Tsurf=-10.0)
+f = heateq_with_phase_change(grid, T_ub)
+jac_prototype = Tridiagonal(ones(length(H0)-1), ones(length(H0)), ones(length(H0)-1))
+tspan = (0.0, 48*3600.0)
+prob = ODEProblem(ODEFunction(f; jac_prototype), H0, tspan, p)
+sol = @time solve(prob, SSPRK43(), abstol=1e-2, reltol=1e-6)
+
+H_sol = reduce(hcat, sol.u)
+T_sol = enthalpyinv.(H_sol, Ref(p))
+Plots.plot(T_sol[1:5:50,:]', leg=nothing)
+
+function loss(p)
+    sensealg = sensealg=InterpolatingAdjoint(autojacvec=EnzymeVJP(), checkpointing=true)
+    sol_p = solve(prob, SSPRK43(); p, sensealg)
+    return mean(sol_p.u[end])/L
+end
+
+grad = @time Zygote.gradient(loss, p)