diff --git a/Project.toml b/Project.toml
index 6fa3f1a..f4d0ae4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,7 +29,7 @@ ProgressMeter = "1"
 PyCall = "1"
 Reexport = "1"
 Revise = "3"
-Sleipnir = "0.2"
+Sleipnir = "0.3"
 Tullio = "0.3"
 
 [extras]
diff --git a/src/models/iceflow/SIA2D/SIA2D.jl b/src/models/iceflow/SIA2D/SIA2D.jl
index 25cb50a..18e739a 100644
--- a/src/models/iceflow/SIA2D/SIA2D.jl
+++ b/src/models/iceflow/SIA2D/SIA2D.jl
@@ -9,6 +9,7 @@ include("SIA2D_utils.jl")
 
 mutable struct SIA2Dmodel{F <: AbstractFloat, I <: Integer} <: SIAmodel
     A::Union{Ref{F}, Nothing}
+    n::Union{Ref{F}, Nothing}
     H::Union{Matrix{F}, Nothing}
     H̄::Union{Matrix{F}, Nothing}
     S::Union{Matrix{F}, Nothing}
@@ -38,6 +39,7 @@ end
 
 function SIA2Dmodel(params::Sleipnir.Parameters;
                     A::Union{Ref{F}, Nothing} = nothing,
+                    n::Union{Ref{F}, Nothing} = nothing,
                     H::Union{Matrix{F}, Nothing} = nothing,
                     H̄::Union{Matrix{F}, Nothing} = nothing,
                     S::Union{Matrix{F}, Nothing} = nothing,
@@ -66,7 +68,7 @@ function SIA2Dmodel(params::Sleipnir.Parameters;
     
     ft = params.simulation.float_type
     it = params.simulation.int_type
-    SIA2D_model = SIA2Dmodel{ft,it}(A, H, H̄, S, dSdx, dSdy, D, Dx, Dy, dSdx_edges, dSdy_edges,
+    SIA2D_model = SIA2Dmodel{ft,it}(A, n, H, H̄, S, dSdx, dSdy, D, Dx, Dy, dSdx_edges, dSdy_edges,
                             ∇S, ∇Sx, ∇Sy, Fx, Fy, Fxx, Fyy, V, Vx, Vy, Γ, MB, MB_mask, MB_total, glacier_idx)
 
     return SIA2D_model
@@ -95,7 +97,8 @@ function initialize_iceflow_model!(iceflow_model::IF,
                                    ) where {IF <: IceflowModel, I <: Int}
     nx, ny = glacier.nx, glacier.ny
     F = params.simulation.float_type
-    iceflow_model.A = Ref{F}(params.physical.A)
+    iceflow_model.A = Ref{F}(glacier.A)
+    iceflow_model.n = Ref{F}(glacier.n)
     iceflow_model.H = deepcopy(glacier.H₀)::Matrix{F}
     iceflow_model.H̄ = zeros(F,nx-1,ny-1)
     iceflow_model.S = deepcopy(glacier.S)::Matrix{F}
diff --git a/src/models/iceflow/SIA2D/SIA2D_utils.jl b/src/models/iceflow/SIA2D/SIA2D_utils.jl
index 02c69d7..4cdb251 100644
--- a/src/models/iceflow/SIA2D/SIA2D_utils.jl
+++ b/src/models/iceflow/SIA2D/SIA2D_utils.jl
@@ -14,6 +14,7 @@ function SIA2D!(dH::Matrix{F}, H::Matrix{F}, simulation::SIM, t::F) where {F <:
         int_type = simulation.parameters.simulation.int_type
     H̄::Matrix{F} = SIA2D_model.H̄
     A::Ref{F} = SIA2D_model.A
+    n::Ref{F} = SIA2D_model.n
     B::Matrix{F} = glacier.B
     S::Matrix{F} = SIA2D_model.S
     dSdx::Matrix{F} = SIA2D_model.dSdx
@@ -33,7 +34,6 @@ function SIA2D!(dH::Matrix{F}, H::Matrix{F}, simulation::SIM, t::F) where {F <:
     Δx::F = glacier.Δx
     Δy::F = glacier.Δy
     Γ::Ref{F} = SIA2D_model.Γ
-    n::int_type = simulation.parameters.physical.n
     ρ::F = simulation.parameters.physical.ρ
     g::F = simulation.parameters.physical.g
 
@@ -48,12 +48,12 @@ function SIA2D!(dH::Matrix{F}, H::Matrix{F}, simulation::SIM, t::F) where {F <:
     diff_y!(dSdy, S, Δy) 
     avg_y!(∇Sx, dSdx)
     avg_x!(∇Sy, dSdy)
-    ∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n - 1)/2) 
+    ∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n[] - 1)/2) 
 
     # @infiltrate
     avg!(H̄, H)
-    Γ[] = 2.0 * A[] * (ρ * g)^n / (n+2) # 1 / m^3 s 
-    D .= Γ[] .* H̄.^(n + 2) .* ∇S
+    Γ[] = 2.0 * A[] * (ρ * g)^n[] / (n[]+2) # 1 / m^3 s 
+    D .= Γ[] .* H̄.^(n[] + 2) .* ∇S
 
     # Compute flux components
     @views diff_x!(dSdx_edges, S[:,2:end - 1], Δx)
@@ -94,6 +94,7 @@ function SIA2D(H, SIA2Dmodel)
     Δx = SIA2Dmodel.Δx
     Δy = SIA2Dmodel.Δy
     A = SIA2Dmodel.A
+    n = SIA2Dmodel.n
 
     # Update glacier surface altimetry
     S = B .+ H
@@ -174,7 +175,7 @@ function avg_surface_V(ifm::IF, temp::F, sim, θ=nothing, UA_f=nothing; testmode
 end
 
 """
-    surface_V(H, B, Δx, Δy, temp, sim, A_noise, θ=[], UA_f=[])
+    surface_V!(H, B, Δx, Δy, temp, sim, A_noise, θ=[], UA_f=[])
 
 Computes the ice surface velocity for a given glacier state
 """
@@ -196,9 +197,9 @@ function surface_V!(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SI
     Dx::Matrix{ft} = iceflow_model.Dx
     Dy::Matrix{ft} = iceflow_model.Dy
     A::Ref{ft} = iceflow_model.A
+    n::Ref{ft} = iceflow_model.n
     Δx::ft = glacier.Δx
     Δy::ft = glacier.Δy
-    n::it = params.physical.n
     ρ::ft = params.physical.ρ
     g::ft = params.physical.g
     
@@ -211,11 +212,11 @@ function surface_V!(H::Matrix{F}, simulation::SIM) where {F <: AbstractFloat, SI
     diff_y!(dSdy, S, Δy) 
     avg_y!(∇Sx, dSdx)
     avg_x!(∇Sy, dSdy)
-    ∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n - 1)/2) 
+    ∇S .= (∇Sx.^2 .+ ∇Sy.^2).^((n[] - 1)/2) 
 
     avg!(H̄, H)
-    Γꜛ[] = 2.0 * A[] * (ρ * g)^n / (n+1) # surface stress (not average)  # 1 / m^3 s 
-    D .= Γꜛ[] .* H̄.^(n + 1) .* ∇S
+    Γꜛ[] = 2.0 * A[] * (ρ * g)^n[] / (n[]+1) # surface stress (not average)  # 1 / m^3 s 
+    D .= Γꜛ[] .* H̄.^(n[] + 1) .* ∇S
     
     # Compute averaged surface velocities
     Vx = .-D .* ∇Sx
diff --git a/test/PDE_UDE_solve.jl b/test/PDE_UDE_solve.jl
index a35b95f..a07e0e3 100644
--- a/test/PDE_UDE_solve.jl
+++ b/test/PDE_UDE_solve.jl
@@ -1,6 +1,7 @@
 
 
-function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::Bool=true) where {F <: AbstractFloat}
+function pde_solve_test(; rtol::F, atol::F, save_refs::Bool=false, MB::Bool=false, fast::Bool=true) where {F <: AbstractFloat}
+
     println("PDE solving with MB = $MB")
     working_dir = joinpath(homedir(), "OGGM/Huginn_tests")
     
@@ -9,9 +10,10 @@ function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::
                                               workers=2,
                                               ice_thickness_source = "Farinotti19"),
                         simulation = SimulationParameters(use_MB=MB,
-                                                        velocities=false,
-                                                        tspan=(2010.0, 2015.0),
-                                                        working_dir = Huginn.root_dir)
+                                                          velocities=false,
+                                                          tspan=(2010.0, 2015.0),
+                                                          working_dir = Huginn.root_dir),
+                        solver = SolverParameters(reltol=1e-12)
                         ) 
 
     ## Retrieving gdirs and climate for the following glaciers
@@ -24,13 +26,6 @@ function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::
         "RGI60-07.00274", "RGI60-07.01323",  "RGI60-01.17316"]
     end
 
-    # Load reference values for the simulation
-    if MB
-        PDE_refs = load(joinpath(Huginn.root_dir, "test/data/PDE/PDE_refs_MB.jld2"))["results"]
-    else
-        PDE_refs = load(joinpath(Huginn.root_dir, "test/data/PDE/PDE_refs_noMB.jld2"))["results"]
-    end
-
     model = Model(iceflow = SIA2Dmodel(params))
 
     # We retrieve some glaciers for the simulation
@@ -51,6 +46,13 @@ function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::
         end
     end
     
+    # Load reference values for the simulation
+    if MB
+        PDE_refs = load(joinpath(Huginn.root_dir, "test/data/PDE/PDE_refs_MB.jld2"))["results"]
+    else
+        PDE_refs = load(joinpath(Huginn.root_dir, "test/data/PDE/PDE_refs_noMB.jld2"))["results"]
+    end
+
     # # Run one epoch of the UDE training
     # θ = zeros(10) # dummy data for the NN
     # UA_f = zeros(10)
@@ -59,6 +61,7 @@ function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::
     # H_V_preds = @time predict_iceflow(θ, UA_f, gdirs, context_batches, UDE_settings, mb_model) # Array{(H_pred, V̄x_pred, V̄y_pred, rgi_id)}
 
     let results=prediction.results
+
     for result in results
         let result=result, test_ref=nothing, UDE_pred=nothing
         for PDE_ref in PDE_refs
@@ -79,19 +82,22 @@ function pde_solve_test(atol::F, save_refs::Bool = false; MB::Bool=false, fast::
             mkdir(test_plot_path)
         end
         MB ? vtol = 30.0*atol : vtol = 12.0*atol # a little extra tolerance for UDE surface velocities
+
         ### PDE ###
-        plot_test_error(result, test_ref, "H", result.rgi_id, atol, MB)
+        plot_test_error(result, test_ref, "H",  result.rgi_id, atol, MB)
         plot_test_error(result, test_ref, "Vx", result.rgi_id, vtol, MB)
         plot_test_error(result, test_ref, "Vy", result.rgi_id, vtol, MB)
+        
         ### UDE ###
         # Huginn.plot_test_error(UDE_pred, test_ref, "H", rgi_id, atol, MB)
         # Huginn.plot_test_error(UDE_pred, test_ref, "Vx", rgi_id, vtol, MB)
         # Huginn.plot_test_error(UDE_pred, test_ref, "Vy", rgi_id, vtol, MB)
 
         # Test that the PDE simulations are correct
-        @test all(isapprox.(result.H[end], test_ref.H[end], atol=atol))
-        @test all(isapprox.(result.Vx, test_ref.Vx, atol=vtol))
-        @test all(isapprox.(result.Vy, test_ref.Vy, atol=vtol))
+        @test all(isapprox.(result.H[end], test_ref.H[end], rtol=rtol, atol=atol))
+        @test all(isapprox.(result.Vx, test_ref.Vx, rtol=rtol, atol=vtol))
+        @test all(isapprox.(result.Vy, test_ref.Vy, rtol=rtol, atol=vtol))
+        
         # Test that the UDE simulations are correct
         # @test all(isapprox.(UDE_pred[1], test_ref["H"], atol=atol))
         # @test all(isapprox.(UDE_pred[2], test_ref["Vx"], atol=vtol)) 
diff --git a/test/data/PDE/PDE_refs_noMB.jld2 b/test/data/PDE/PDE_refs_noMB.jld2
index 7ce09c1..69a1e8b 100644
Binary files a/test/data/PDE/PDE_refs_noMB.jld2 and b/test/data/PDE/PDE_refs_noMB.jld2 differ
diff --git a/test/halfar.jl b/test/halfar.jl
index 10244a2..f073870 100644
--- a/test/halfar.jl
+++ b/test/halfar.jl
@@ -17,7 +17,7 @@ Arguments
     - `distance_to_border`: Minimum distance to the border used to evaluate test. Points close to the border are not considered. 
     - `save_plot`: Save plot with comparision of prediction and true solution.
 """
-function unit_halfar_test(; A, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol=0.02, atol=1.0, distance_to_border=3, save_plot=false)
+function unit_halfar_test(; A, n, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol=0.02, atol=1.0, distance_to_border=3, save_plot=false)
 
     # Get parameters for a simulation 
     parameters = Parameters(simulation=SimulationParameters(tspan=(t₀, t₁),
@@ -26,26 +26,24 @@ function unit_halfar_test(; A, t₀, t₁, Δx, Δy, nx, ny, h₀, r₀, rtol=0.
                                                             multiprocessing=true,
                                                             workers=1,
                                                             working_dir=Huginn.root_dir),
-                            physical=PhysicalParameters(A=A),
+                            physical=PhysicalParameters(),
                             solver=SolverParameters(reltol=1e-12))
 
     # Bed (it has to be flat for the Halfar solution)
     B = zeros((nx,ny))
 
-    model = Model(iceflow = SIA2Dmodel(parameters)) #,
-                #   mass_balance = nothing, # TImodel1(parameters; DDF=8.0/1000.0, acc_factor=1.0/1000.0),
-                #   machine_learning = nothing) # NN(parameters))
+    model = Model(iceflow = SIA2Dmodel(parameters)) 
 
     # Initial condition of the glacier
     R₀ = [sqrt((Δx * (i - nx/2))^2 + (Δy * (j - ny/2))^2) for i in 1:nx, j in 1:ny]
-    H₀ = halfar_solution(R₀, t₀, h₀, r₀, parameters.physical)
+    H₀ = halfar_solution(R₀, t₀, h₀, r₀, A, n, parameters.physical)
     S = B + H₀
     # Final expected solution
-    H₁ = halfar_solution(R₀, t₁, h₀, r₀, parameters.physical)
+    H₁ = halfar_solution(R₀, t₁, h₀, r₀, A, n, parameters.physical)
 
     # Define glacier object
-    glacier = Glacier(rgi_id = "toy", H₀ = H₀, S = S, B = B, 
-                      Δx=Δx, Δy=Δy, nx=nx, ny=ny)
+    glacier = Glacier2D(rgi_id = "toy", H₀ = H₀, S = S, B = B, A = A, n=n, 
+                        Δx=Δx, Δy=Δy, nx=nx, ny=ny)
     glaciers = Vector{Sleipnir.AbstractGlacier}([glacier])
 
     prediction = Prediction(model, glaciers, parameters)
@@ -93,13 +91,13 @@ Arguments
     - `atol`: Absolute tolerance for ice thickness H 
 """
 function halfar_test(; rtol, atol)
-    unit_halfar_test(A=4e-17, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol)
-    unit_halfar_test(A=8e-17, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol)
-    unit_halfar_test(A=4e-17, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
-    unit_halfar_test(A=8e-17, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
-    unit_halfar_test(A=4e-17, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
-    unit_halfar_test(A=8e-17, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
-    unit_halfar_test(A=4e-17, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol)
-    unit_halfar_test(A=8e-17, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol)
-    unit_halfar_test(A=4e-17, t₀=5.0, t₁=10.0, Δx=10.0, Δy=10.0, nx=500, ny=500, h₀=300, r₀=1000, rtol=rtol, atol=atol)
+    unit_halfar_test(A=4e-17, n=3.0, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol)
+    unit_halfar_test(A=8e-17, n=3.0, t₀=5.0, t₁=10.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=1000, rtol=rtol, atol=atol)
+    unit_halfar_test(A=4e-17, n=3.0, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
+    unit_halfar_test(A=8e-17, n=3.0, t₀=5.0, t₁=40.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
+    unit_halfar_test(A=4e-17, n=3.0, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
+    unit_halfar_test(A=8e-17, n=3.0, t₀=5.0, t₁=100.0, Δx=50.0, Δy=50.0, nx=100, ny=100, h₀=500, r₀=600, rtol=rtol, atol=atol)
+    unit_halfar_test(A=4e-17, n=3.0, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol)
+    unit_halfar_test(A=8e-17, n=3.0, t₀=5.0, t₁=40.0, Δx=80.0, Δy=80.0, nx=100, ny=100, h₀=300, r₀=1000, rtol=rtol, atol=atol)
+    unit_halfar_test(A=4e-17, n=3.0, t₀=5.0, t₁=10.0, Δx=10.0, Δy=10.0, nx=500, ny=500, h₀=300, r₀=1000, rtol=rtol, atol=atol)
 end
\ No newline at end of file
diff --git a/test/mass_conservation.jl b/test/mass_conservation.jl
index fbf885e..93affe0 100644
--- a/test/mass_conservation.jl
+++ b/test/mass_conservation.jl
@@ -21,7 +21,7 @@ Arguments
     - `rtol`: Relative tolerance
     - `save_plot`: Optional bool to save plot during simulation
 """
-function unit_mass_test(; H₀, B, A, t_sim, Δx, Δy, rtol=0.02, save_plot=false)
+function unit_mass_test(; H₀, B, A, n, t_sim, Δx, Δy, rtol=0.02, save_plot=false)
 
     # Get parameters for a simulation 
     parameters = Parameters(simulation=SimulationParameters(tspan=(0.0, t_sim),
@@ -29,7 +29,7 @@ function unit_mass_test(; H₀, B, A, t_sim, Δx, Δy, rtol=0.02, save_plot=fals
                                                             use_iceflow=true,
                                                             multiprocessing=true,
                                                             workers=1),
-                            physical=PhysicalParameters(A=A),
+                            physical=PhysicalParameters(),
                             solver=SolverParameters(reltol=1e-12))
 
     model = Model(iceflow = SIA2Dmodel(parameters))
@@ -39,8 +39,8 @@ function unit_mass_test(; H₀, B, A, t_sim, Δx, Δy, rtol=0.02, save_plot=fals
     nx, ny = size(H₀)
 
     # Define glacier object
-    glacier = Glacier(rgi_id = "toy", H₀ = H₀, S = S, B = B, 
-                      Δx=Δx, Δy=Δy, nx=nx, ny=ny)
+    glacier = Glacier2D(rgi_id = "toy", H₀ = H₀, S = S, B = B, A = A, n = n, 
+                        Δx=Δx, Δy=Δy, nx=nx, ny=ny)
     glaciers = Vector{Sleipnir.AbstractGlacier}([glacier])
 
     prediction = Prediction(model, glaciers, parameters)
@@ -99,7 +99,7 @@ function unit_mass_flatbed_test(; rtol)
                     H₀ = zeros((nx,ny))
                     @views H₀[floor(Int,nx/3):floor(Int,2nx/3), floor(Int,ny/3):floor(Int,2ny/3)] .= 400
                 end
-                unit_mass_test(; H₀=H₀, B=B, A=A, t_sim=10.0, Δx=50.0, Δy=50.0, rtol=rtol, save_plot=false)
+                unit_mass_test(; H₀=H₀, B=B, A=A, n=3.0, t_sim=10.0, Δx=50.0, Δy=50.0, rtol=rtol, save_plot=false)
             end
         end
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 89dc978..9d3802d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,10 +23,9 @@ ENV["GKSwstype"]="nul"
 
 @testset "Solver parameters construction with default variables" params_constructor_default()
 
-@testset "Halfar Solution" halfar_test(rtol=0.02, atol=1.0)
+@testset "Halfar Solution" halfar_test(; rtol=0.02, atol=1.0)
 
-atol = 0.01
-@testset "PDE solving integration tests" pde_solve_test(atol; MB=false, fast=true)
+@testset "PDE solving integration tests" pde_solve_test(; rtol=0.01, atol=0.01, save_refs=false, MB=false, fast=true)
 
-@testset "Conservation of Mass - Flat Bed" unit_mass_flatbed_test(rtol=1.0e-7)
+@testset "Conservation of Mass - Flat Bed" unit_mass_flatbed_test(; rtol=1.0e-7)
 
diff --git a/test/utils_test.jl b/test/utils_test.jl
index bf4f0a7..acb9e3a 100644
--- a/test/utils_test.jl
+++ b/test/utils_test.jl
@@ -27,13 +27,11 @@ Arguments:
     - t: time
     - ν = (A, H₀, R₀) 
 """
-function halfar_solution(R, t, h₀, r₀, physical_parameters::PhysicalParameters)
+function halfar_solution(R, t, h₀, r₀, A, n, physical_parameters::PhysicalParameters)
 
     # parameters of Halfar solutions
     ρ = physical_parameters.ρ
     g = physical_parameters.g
-    n = physical_parameters.n
-    A = physical_parameters.A
 
     Γ = 2 * A * (ρ * g)^n / (n+2)
     # Characteristic time