diff --git a/examples/example_config_files/ivc_inactivelayer.yaml b/examples/example_config_files/ivc_inactivelayer.yaml
index b99ee6378..3bc9aff87 100644
--- a/examples/example_config_files/ivc_inactivelayer.yaml
+++ b/examples/example_config_files/ivc_inactivelayer.yaml
@@ -20,6 +20,7 @@ grid:
       boundaries:
         left: inf
         right: inf
+  spacing_surface_refinement: [1e-4, 1e-4, 1e-4] #m
 medium: vacuum
 detectors:
 - semiconductor:
diff --git a/examples/example_config_files/public_CGD_config_cyl_grid.yaml b/examples/example_config_files/public_CGD_config_cyl_grid.yaml
index 999663a1b..c03767cac 100644
--- a/examples/example_config_files/public_CGD_config_cyl_grid.yaml
+++ b/examples/example_config_files/public_CGD_config_cyl_grid.yaml
@@ -10,7 +10,7 @@ grid:
     r:
       to: 9
       boundaries: inf
-    y:
+    phi:
       from: 0
       to: 0
       boundaries: periodic
diff --git a/examples/example_config_files/test_detector.yaml b/examples/example_config_files/test_detector.yaml
index b50c18d08..67d2783e2 100644
--- a/examples/example_config_files/test_detector.yaml
+++ b/examples/example_config_files/test_detector.yaml
@@ -10,7 +10,7 @@ grid:
     r:
       to: 9
       boundaries: inf
-    y:
+    phi:
       from: 0
       to: 180
       boundaries: periodic
diff --git a/src/PotentialCalculation/Refinement.jl b/src/PotentialCalculation/Refinement.jl
index fab5dac6b..9ef1a94de 100644
--- a/src/PotentialCalculation/Refinement.jl
+++ b/src/PotentialCalculation/Refinement.jl
@@ -174,3 +174,76 @@ end
 
 _extend_refinement_limits(rl::Real) = (rl, rl, rl )
 _extend_refinement_limits(rl::Tuple{<:Real,<:Real,<:Real}) = rl
+
+@inline function has_surface_points(slice::AbstractArray{PointType})::Bool
+    return any(is_in_inactive_layer, slice)
+end
+
+function _refine_axis_surface( ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T}}, surface_intervals::AbstractVector{Bool}, min_spacing::T;
+    extra_before::Int = 5,  # intervals to refine before first surface interval
+    extra_after::Int = 5    # intervals to refine after last surface interval
+) where {T}
+
+    old_ticks = ax.ticks
+    n_int = length(surface_intervals)
+
+    # Find first and last surface intervals
+    first_surface = findfirst(surface_intervals)
+    last_surface  = findlast(surface_intervals)
+
+    (isnothing(first_surface) || isnothing(last_surface)) && return ax
+
+    # Create flags for all intervals to refine
+    refine_flags = copy(surface_intervals)
+
+    # Add extra intervals before first surface interval
+    start_idx = max(first_surface - extra_before, 1)
+    refine_flags[start_idx:first_surface-1] .= true
+
+    # Add extra intervals after last surface interval
+    end_idx = min(last_surface + extra_after, n_int)
+    refine_flags[last_surface+1:end_idx] .= true
+
+    # Merge consecutive intervals to refine
+    merged = Vector{UnitRange{Int}}()
+    i = 1
+    while i <= n_int
+        if refine_flags[i]
+            j = i
+            while j < n_int && refine_flags[j+1]
+                j += 1
+            end
+            push!(merged, i:j)
+            i = j + 1
+        else
+            i += 1
+        end
+    end
+
+    # Compute number of points to add per interval
+    ns = zeros(Int, n_int)
+    for r in merged
+        for i in r
+            Δ = old_ticks[i+1] - old_ticks[i]
+            if Δ > min_spacing
+                ns[i] = ceil(Int, Δ/min_spacing) - 1
+            end
+        end
+    end
+
+    # Subdivide intervals
+    sub_widths = [(old_ticks[i+1] - old_ticks[i]) / (ns[i]+1) for i in 1:n_int]
+    ticks = Vector{T}(undef, length(old_ticks) + sum(ns))
+    i_tick = 1
+    for j in 1:n_int
+        ticks[i_tick] = old_ticks[j]
+        for k in 1:ns[j]
+            i_tick += 1
+            ticks[i_tick] = old_ticks[j] + k*sub_widths[j]
+        end
+        i_tick += 1
+    end
+    ticks[end] = old_ticks[end]
+
+    return typeof(ax)(ax.interval, ticks)
+end
diff --git a/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl b/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl
index e6bc1476c..fcb906309 100644
--- a/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl
+++ b/src/PotentialCalculation/SuccessiveOverRelaxation/CPU_innerloop.jl
@@ -10,14 +10,13 @@
     pwΔmp2r, pwΔmp2l, 
     pwΔmp3r, pwΔmp3l,
 ) where {T, S}
-    
 
     # pww1r        = pcs.geom_weights[3][1, in1]
     pww1r        = geom_weights_3[1, in1]
     pww1l        = geom_weights_3[2, in1]
     pwΔmp1       = geom_weights_3[3, in1]
     Δ1_ext_inv_l = geom_weights_3[4, in1]
-    Δ1_ext_inv_r = geom_weights_3[4, in1 + 1]
+    Δ1_ext_inv_r = in1 < size(geom_weights_3, 2) ? geom_weights_3[4, in1 + 1] : geom_weights_3[4, in1]
 
     # ϵ_ijk: i: 1.RB-Dim. | j: 2.RB-Dim. | k: 3.RB-Dim.
     ϵ_lll, ϵ_llr, ϵ_lrl, ϵ_lrr, ϵ_rll, ϵ_rlr, ϵ_rrl, ϵ_rrr = get_ϵ_of_oktant(
@@ -168,18 +167,25 @@ end
     ϵ_r::AbstractArray{T, 3}, ::Type{Cylindrical},
     i1, in1, i2, in2, i3, in3
 ) where {T}
+    # Cylindrical: inner loop over z (3rd dimension)
+    i1_safe  = clamp(i1, 1, size(ϵ_r, 3))
+    in1_safe = clamp(in1, 1, size(ϵ_r, 3))
+    i2_safe  = clamp(i2, 1, size(ϵ_r, 2))
+    in2_safe = clamp(in2, 1, size(ϵ_r, 2))
+    i3_safe  = clamp(i3, 1, size(ϵ_r, 1))
+    in3_safe = clamp(in3, 1, size(ϵ_r, 1))
     # ϵ_r is not transformed into an red-black-4D-array.
     # The inner loop (over i1) is along the z-Dimension (Cylindrical Case), 
     # which is the 3rd dimension for Cylindrical coordinates: (r, φ, z)
     return @inbounds begin
-        ϵ_r[ in3, in2, in1 ],
-        ϵ_r[  i3, in2, in1 ],
-        ϵ_r[ in3,  i2, in1 ],
-        ϵ_r[  i3,  i2, in1 ],
-        ϵ_r[ in3, in2,  i1 ],
-        ϵ_r[  i3, in2,  i1 ],
-        ϵ_r[ in3,  i2,  i1 ],
-        ϵ_r[  i3,  i2,  i1 ]
+        ϵ_r[in3_safe, in2_safe, in1_safe],
+        ϵ_r[i3_safe,  in2_safe, in1_safe],
+        ϵ_r[in3_safe,  i2_safe, in1_safe],
+        ϵ_r[i3_safe,   i2_safe, in1_safe],
+        ϵ_r[in3_safe, in2_safe,  i1_safe],
+        ϵ_r[i3_safe,  in2_safe,  i1_safe],
+        ϵ_r[in3_safe,  i2_safe,  i1_safe],
+        ϵ_r[i3_safe,   i2_safe,  i1_safe]
     end
 end
 
@@ -187,17 +193,24 @@ end
     ϵ_r::AbstractArray{T, 3}, ::Type{Cartesian},
     i1, in1, i2, in2, i3, in3
 ) where {T}
+    # Cartesian: inner loop over x (1st dimension)
+    i1_safe  = clamp(i1, 1, size(ϵ_r, 1))
+    in1_safe = clamp(in1, 1, size(ϵ_r, 1))
+    i2_safe  = clamp(i2, 1, size(ϵ_r, 2))
+    in2_safe = clamp(in2, 1, size(ϵ_r, 2))
+    i3_safe  = clamp(i3, 1, size(ϵ_r, 3))
+    in3_safe = clamp(in3, 1, size(ϵ_r, 3))
     # The inner loop (over i1) is along the x-Dimension (Cartesian Case), 
     # which is the 1rd dimension for Cartesian coordinates: (x, y, z)
     return @inbounds begin
-        ϵ_r[ in1, in2, in3 ], 
-        ϵ_r[ in1, in2,  i3 ],
-        ϵ_r[ in1,  i2, in3 ], 
-        ϵ_r[ in1,  i2,  i3 ],
-        ϵ_r[  i1, in2, in3 ],
-        ϵ_r[  i1, in2,  i3 ],
-        ϵ_r[  i1,  i2, in3 ],
-        ϵ_r[  i1,  i2,  i3 ]
+        ϵ_r[in1_safe, in2_safe, in3_safe],
+        ϵ_r[in1_safe, in2_safe,  i3_safe],
+        ϵ_r[in1_safe,  i2_safe, in3_safe],
+        ϵ_r[in1_safe,  i2_safe,  i3_safe],
+        ϵ_r[ i1_safe, in2_safe, in3_safe],
+        ϵ_r[ i1_safe, in2_safe,  i3_safe],
+        ϵ_r[ i1_safe,  i2_safe, in3_safe],
+        ϵ_r[ i1_safe,  i2_safe,  i3_safe]
     end
 end
 
@@ -242,4 +255,4 @@ end
     np::NTuple{6, T}, ::Type{Cartesian}, i::Int
 ) where {T}
     return np
-end
\ No newline at end of file
+end
diff --git a/src/Simulation/Simulation.jl b/src/Simulation/Simulation.jl
index a8cf25af9..161bc3906 100644
--- a/src/Simulation/Simulation.jl
+++ b/src/Simulation/Simulation.jl
@@ -178,7 +178,7 @@ function Simulation{T}(dict::AbstractDict)::Simulation{T} where {T <: SSDFloat}
     sim.medium = material_properties[materials[haskey(dict, "medium") ? dict["medium"] : "vacuum"]]
     sim.detector = SolidStateDetector{T}(dict, sim.input_units) 
     sim.world = if haskey(dict, "grid") && isa(dict["grid"], AbstractDict) && haskey(dict["grid"], "axes")
-            World(T, dict["grid"], sim.input_units)
+        World(T, dict["grid"], sim.input_units)
         else let det = sim.detector 
             world_limits = get_world_limits_from_objects(CS, det)
             World(CS, world_limits)
@@ -759,6 +759,7 @@ Takes the current state of `sim.electric_potential` and refines it with respect
 SolidStateDetectors.refine!(sim, ElectricPotential, max_diffs = (100, 100, 100), minimum_distances = (0.01, 0.02, 0.01))
 ```
 """
+
 function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
                     max_diffs::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0)),
                     minimum_distances::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0));
@@ -768,9 +769,10 @@ function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
     sim.electric_potential = refine_scalar_potential(sim.electric_potential, T.(max_diffs), T.(minimum_distances))
 
     if update_other_fields
-        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data,
-                                        not_only_paint_contacts = not_only_paint_contacts, paint_contacts = paint_contacts)
+        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data;
+                                        not_only_paint_contacts, paint_contacts)
 
+        
         sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
         sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
         sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
@@ -812,6 +814,123 @@ function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int
     nothing
 end
 
+"""
+        refine_surface!(sim::Simulation{T};
+                    min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4)),
+                    not_only_paint_contacts::Bool = true,
+                    paint_contacts::Bool = true,
+                    update_other_fields::Bool = true) where {T <: SSDFloat}
+
+Refine the simulation's electric potential grid **only in regions containing surface points**, keeping other regions unchanged.
+
+This function:
+
+1. Identifies which intervals along each axis contain surface points.
+2. Refines these intervals by adding extra grid points while respecting `min_spacing`.
+3. Interpolates the existing electric potential onto the refined grid.
+4. Updates the simulation's `electric_potential` and optionally dependent fields (`imp_scale`, `q_eff_imp`, `q_eff_fix`, `ϵ_r`, `point_types`).
+
+Special behaviour:
+
+- In cylindrical coordinates, the φ-axis (2nd axis) is left unchanged except for closing the axis for interpolation purposes.
+- Handles both 2D and 3D grids.
+
+# Arguments
+- `sim::Simulation{T}`: The simulation object containing the electric potential to refine.
+- `min_spacing::NTuple{3,T}`: Minimum allowed spacing in the surface after refinement.
+"""
+
+function refine_surface!(sim::Simulation{T,CS}, min_spacing::NTuple{3,T} = (T(1e-4), T(1e-4), T(1e-4));
+                         not_only_paint_contacts::Bool = true,
+                         paint_contacts::Bool = true,
+                         update_other_fields::Bool = true) where {T <: SSDFloat, CS <: AbstractCoordinateSystem}
+    
+    old_grid = sim.electric_potential.grid
+
+    # Enforce minimum and maximum spacing for surface refinement (10 µm < min_spacing < 1 mm)
+    for d in 1:3
+        if length(old_grid.axes[d].ticks) != 1
+            if min_spacing[d] < T(1e-5) || min_spacing[d] > T(1e-3)
+                @warn "min_spacing[$d] = $(min_spacing[d]) m is out of bounds (1e-5 m < min_spacing < 1e-3 m). Higher values don't require surface refinement"
+                if min_spacing[d] < T(1e-5)
+                    @warn "Using min_spacing[$d] = 1e-5 m"
+                else
+                    @warn "Using min_spacing[$d] = 1e-3 m"  
+                end
+            end
+        end
+    end
+    
+    min_spacing = ntuple(d -> length(old_grid.axes[d].ticks) == 1 ? T(0.0) : clamp(T(min_spacing[d]), T(1e-5), T(1e-3)), 3)
+
+    # Determine which intervals have surface points
+    surface_intervals = ntuple(d -> begin
+        n_int = length(old_grid.axes[d].ticks) - 1
+        map(i -> begin
+            slice = if d == 1
+                view(sim.point_types.data, i, :, :)
+            elseif d == 2
+                view(sim.point_types.data, :, i, :)
+            else
+                view(sim.point_types.data, :, :, i)
+            end
+            has_surface_points(slice)
+        end, 1:n_int)
+    end, 3)
+
+    # Refine axes: skip phi-axis if Cylindrical
+    is_cyl = CS === Cylindrical
+    new_axes = ntuple(i -> begin
+                          if is_cyl && i == 2
+                              # Keep φ-axis type but close it for interpolation
+                              _get_closed_axis(old_grid.axes[i])
+                          else
+                              _refine_axis_surface( old_grid.axes[i], surface_intervals[i], min_spacing[i],)
+                          end
+                      end, 3)
+
+    new_grid = Grid{T, 3, CS}((new_axes[1], new_axes[2], new_axes[3]))
+    
+    # Interpolate potential onto new grid
+    closed_pot = _get_closed_potential(sim.electric_potential)
+    new_data = Array{T,3}(undef, size(new_grid))
+    only2d = is_cyl && size(closed_pot.data, 2) == 1
+    int = interpolate_closed_potential(closed_pot, Val(only2d))
+    
+    if only2d
+        # 2D grid
+        for i3 in axes(new_data,3), i1 in axes(new_data,1)
+            new_data[i1,1,i3] = int(new_grid.axes[1].ticks[i1], new_grid.axes[3].ticks[i3])
+        end
+    else
+        # 3D grid
+        for i3 in axes(new_data,3), i2 in axes(new_data,2), i1 in axes(new_data,1)
+            new_data[i1,i2,i3] = int(new_grid.axes[1].ticks[i1],
+                                     new_grid.axes[2].ticks[i2],
+                                     new_grid.axes[3].ticks[i3])
+        end
+    end
+    
+    # Update electric potential in simulation
+    sim.electric_potential = _convert_to_original_potential(sim.electric_potential, new_data, new_grid)
+    
+    # Update dependent fields
+    if update_other_fields
+
+        pcs = PotentialCalculationSetup(sim.detector, sim.electric_potential.grid, sim.medium, sim.electric_potential.data;
+                                        not_only_paint_contacts, paint_contacts)
+        sim.imp_scale = ImpurityScale(ImpurityScaleArray(pcs), sim.electric_potential.grid)
+        sim.q_eff_imp = EffectiveChargeDensity(EffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
+        sim.q_eff_fix = EffectiveChargeDensity(FixedEffectiveChargeDensityArray(pcs), sim.electric_potential.grid)
+        sim.ϵ_r = DielectricDistribution(DielectricDistributionArray(pcs), get_extended_midpoints_grid(sim.electric_potential.grid))
+        sim.point_types = PointTypes(PointTypeArray(pcs), sim.electric_potential.grid)
+    end
+
+    return nothing
+end
+
+
+
 """
     compute_min_tick_distance(grid)
 
@@ -899,7 +1018,26 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
         else
             sor_consts = T.(sor_consts)
         end
-      
+        has_surface_model = isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
+
+        if has_surface_model
+            if isnothing(sim.world.spacing_surface_refinement)
+                @warn """Surface model detected but `spacing_surface_refinement` is not defined.
+                         Surface refinement will not be performed."""
+            else
+                if length(refinement_limits) < 3 || last(refinement_limits) > 0.05
+                    @warn """Surface model detected:
+                             - Number of refinement steps: $(length(refinement_limits))
+                             - Last refinement: $(last(refinement_limits))
+                             Surface refinement requires:
+                             - At least 3 refinement steps
+                             - Last refinement ≤ 0.05"""
+                    refinement_limits = [0.2, 0.1, 0.05]
+                    @warn "Falling back to default refinement_limits = $(refinement_limits)"
+                end
+            end    
+        end
+        
         new_min_tick_distance::NTuple{3,T} = begin
             if ismissing(min_tick_distance)
                 compute_min_tick_distance(grid)
@@ -978,6 +1116,9 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                     device_array_type = device_array_type,
                                     use_nthreads = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid), max_nthreads[1], CS) : max_nthreads[1],
                                     sor_consts = sor_consts )
+
+
+            
         else
             apply_initial_state!(sim, potential_type, contact_id, grid; not_only_paint_contacts, paint_contacts, depletion_handling)
             update_till_convergence!( sim, potential_type, contact_id, convergence_limit,
@@ -1003,6 +1144,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                 else
                     abs.(ref_limits .* bias_voltage)
                 end
+                
                 refine!(sim, ElectricPotential, max_diffs, new_min_tick_distance)
                 nt = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid), max_nthreads[iref+1], CS) : max_nthreads[iref+1]
                 verbose && println("Grid size: $(size(sim.electric_potential.data)) - $(onCPU ? "using $(nt) threads now" : "GPU")") 
@@ -1015,6 +1157,23 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
                                                 not_only_paint_contacts = not_only_paint_contacts, 
                                                 paint_contacts = paint_contacts,
                                                 sor_consts = is_last_ref ? T(1) : sor_consts )
+                
+                if has_surface_model && iref==3 && !isnothing(sim.world.spacing_surface_refinement)
+                    mark_bulk_bits!(sim.point_types.data)
+                    mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
+                    mark_inactivelayer_bits!(sim.point_types.data)
+                    verbose && println("Surface Refinement")
+                    # Maximum spacing between (refined) surface ticks (if min_spacing = 1e-4m, only surface intervals wider than 0.1mm get new ticks)
+                    refine_surface!(sim, sim.world.spacing_surface_refinement; update_other_fields=true)
+
+                    update_till_convergence!( sim, potential_type, convergence_limit,
+                                              n_iterations_between_checks = n_iterations_between_checks,
+                                              max_n_iterations = max_n_iterations,
+                                              depletion_handling = depletion_handling,
+                                              device_array_type = device_array_type,
+                                              use_nthreads = guess_nt ? _guess_optimal_number_of_threads_for_SOR(size(sim.electric_potential.grid),
+                                              max_nthreads[1], CS) : max_nthreads[1], sor_consts = T(1) )
+                end
             else
                 max_diffs = abs.(ref_limits)
                 refine!(sim, WeightingPotential, contact_id, max_diffs, new_min_tick_distance)
@@ -1057,7 +1216,7 @@ function _calculate_potential!( sim::Simulation{T, CS}, potential_type::UnionAll
     if depletion_handling && isEP
         mark_undep_bits!(sim.point_types.data, sim.imp_scale.data)
         
-        if isdefined(sim.detector.semiconductor.impurity_density_model, :surface_imp_model)
+        if has_surface_model
             mark_inactivelayer_bits!(sim.point_types.data)
         end
     end
diff --git a/src/World/World.jl b/src/World/World.jl
index 093c15c9d..b382b0d24 100644
--- a/src/World/World.jl
+++ b/src/World/World.jl
@@ -55,18 +55,25 @@ Definition of the finite volume on which a [`Simulation`](@ref) is performed.
 """
 struct World{T <: SSDFloat, N, S} <: AbstractWorld{T, N} 
     intervals::NTuple{N, SSDInterval{T}}
+    spacing_surface_refinement::Union{NTuple{N,T}, Nothing}
 end
 
 function World{T, N, S}(args...) where {T <: SSDFloat, N, S} 
     return World{T, N, S}(args)
 end
 
-
 function World(T, dict::AbstractDict, input_units::NamedTuple)::World
+    
+    spacing_surface_refinement = if haskey(dict, "spacing_surface_refinement")
+        ntuple(i -> _parse_value(T, dict["spacing_surface_refinement"][i], internal_length_unit), 3)
+    else
+        nothing
+    end
+    
     if dict["coordinates"] == "cylindrical"
-        CylindricalWorld(T, dict["axes"], input_units)
+        return CylindricalWorld(T, dict["axes"], input_units, spacing_surface_refinement)
     elseif dict["coordinates"] == "cartesian"
-        CartesianWorld(T, dict["axes"], input_units)
+        return CartesianWorld(T, dict["axes"], input_units, spacing_surface_refinement)
     else
         error("Gridtype must be \"cylindrical\" or \"cartesian\"")
     end
@@ -158,37 +165,37 @@ function get_cartesian_SSDInterval(T, dict::AbstractDict, input_units::NamedTupl
 end
 
 
-function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple)::World
+function CylindricalWorld(T, dict::AbstractDict, input_units::NamedTuple, spacing_surface_refinement)::World
     r_int = get_r_SSDInterval(T, dict["r"], input_units)
     φ_int = get_φ_SSDInterval(T, dict, input_units)
     z_int = get_cartesian_SSDInterval(T, dict["z"], input_units)
-    return World{T, 3, Cylindrical}( r_int, φ_int, z_int )
+    return World{T, 3, Cylindrical}( (r_int, φ_int, z_int), spacing_surface_refinement )
 end
 
 
-function CartesianWorld(T, dict::AbstractDict, input_units::NamedTuple)::World
+function CartesianWorld(T, dict::AbstractDict, input_units::NamedTuple, spacing_surface_refinement)::World
     x_int = get_cartesian_SSDInterval(T, dict["x"], input_units)
     y_int = get_cartesian_SSDInterval(T, dict["y"], input_units)
     z_int = get_cartesian_SSDInterval(T, dict["z"], input_units)
-    return World{T, 3, Cartesian}( x_int, y_int, z_int )
+    return World{T, 3, Cartesian}( (x_int, y_int, z_int), spacing_surface_refinement )
 end
 
-function CartesianWorld(xl::T, xr::T, yl::T, yr::T, zl::T, zr::T)::World where {T <: SSDFloat}
+function CartesianWorld(xl::T, xr::T, yl::T, yr::T, zl::T, zr::T; spacing_surface_refinement=nothing)::World where {T <: SSDFloat}
     Δx::T = (xr - xl) / 10
     Δy::T = (yr - yl) / 10
     Δz::T = (zr - zl) / 10
     x_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(xl - Δx, xr + Δx)
     y_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(yl - Δy, yr + Δy)
     z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - Δz, zr + Δz)
-    return World{T, 3, Cartesian}( x_int, y_int, z_int )
+    return World{T, 3, Cartesian}( (x_int, y_int, z_int), spacing_surface_refinement )
 end
 
-function CylindricalWorld(r_max::T, zl::T, zr::T)::World where {T <: SSDFloat}
+function CylindricalWorld(r_max::T, zl::T, zr::T; spacing_surface_refinement=nothing)::World where {T <: SSDFloat}
     r_int = SSDInterval{T, :closed, :closed, :r0, :infinite}(T(0), abs(r_max * T(1.1)))
     φ_int = SSDInterval{T, :closed, :open, :periodic, :periodic}(T(0), T(2π))
     Δz::T = zr - zl
     z_int = SSDInterval{T, :closed, :closed, :infinite, :infinite}(zl - (Δz / 10), zr + (Δz / 10))
-    return World{T, 3, Cylindrical}( r_int, φ_int, z_int )
+    return World{T, 3, Cylindrical}( (r_int, φ_int, z_int), spacing_surface_refinement )
 end
 
 function World(::Type{Cylindrical}, limits::NTuple{6, T})::World where {T <: SSDFloat}
@@ -205,4 +212,4 @@ end
 function max_tick_distance_default(w::World{T, 3, Cartesian}; n = 50) where {T}
     Δw = width.(w.intervals)
     return max(Δw...) * internal_length_unit / n
-end
\ No newline at end of file
+end
diff --git a/test/test_refine_existing_potential.jl b/test/test_refine_existing_potential.jl
index c65bd0ff6..2dd34ce36 100644
--- a/test/test_refine_existing_potential.jl
+++ b/test/test_refine_existing_potential.jl
@@ -120,3 +120,105 @@ end
     @test length(sim_cart_mm.electric_potential.grid[2]) == 60
     @test length(sim_cart_mm.electric_potential.grid[3]) == 64
 end
+
+@testset "Surface refinement edge cases" begin
+    T = Float64
+    # Surface model but no spacing_surface_refinement
+    begin
+        sim = @test_nowarn Simulation{T}(SSD_examples[:TrueCoaxial])
+        
+        @test_logs (:warn, r"Surface model detected but `spacing_surface_refinement` is not defined") begin
+            timed_calculate_electric_potential!(sim, verbose = false, depletion_handling = true)
+        end
+    end
+    
+    # Wrong refinement limits
+    begin
+        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+        config_dict["grid"]["spacing_surface_refinement"] = [1e-3, 1e-3, 1e-3]
+        sim_cyl = @test_nowarn Simulation{T}(config_dict)
+
+        # Test normal behaviour Cylindrical
+        timed_calculate_electric_potential!(sim_cyl, verbose = false, depletion_handling = true)
+        grid_ax1_cyl = length(sim_cyl.electric_potential.grid[1])
+        grid_ax2_cyl = length(sim_cyl.electric_potential.grid[2])
+        grid_ax3_cyl = length(sim_cyl.electric_potential.grid[3])
+        
+        # --- Case: < 3 refinements ---
+        bad_limits = [0.2, 0.1]
+        @test_logs (:warn, r"Surface model detected") (:warn, r"Falling back to default") match_mode=:any begin
+            timed_calculate_electric_potential!(sim_cyl, refinement_limits=bad_limits, verbose = false, depletion_handling = true)
+        end
+        
+        # --- Case: last refinement > 0.05 ---
+        bad_limits = [0.2, 0.1, 0.09]
+        @test_logs (:warn, r"Surface model detected:") (:warn, r"Falling back to default") match_mode=:any begin
+            timed_calculate_electric_potential!(sim_cyl, refinement_limits=bad_limits, verbose = false, depletion_handling = true)
+        end
+
+        # Test normal behaviour Cartesian
+        config_dict["grid"]["coordinates"] = "cartesian"
+        config_dict["grid"]["axes"]["x"] = Dict(
+            "from" => "-40",
+            "to"   => "40",
+            "boundaries"   => "inf"
+        )
+        config_dict["grid"]["axes"]["y"] = Dict(
+            "from" => "-40",
+            "to"   => "40",
+            "boundaries"   => "inf"
+        )
+        config_dict["grid"]["axes"]["z"] = Dict(
+            "from" => "-10",
+            "to"   => "90",
+            "boundaries"   => "inf"
+        )
+        sim_cart = @test_nowarn Simulation{T}(config_dict)
+        timed_calculate_electric_potential!(sim_cart, verbose = false, depletion_handling = true)
+        grid_ax1_cart = length(sim_cart.electric_potential.grid[1])
+        grid_ax2_cart = length(sim_cart.electric_potential.grid[2])
+	grid_ax3_cart = length(sim_cart.electric_potential.grid[3])
+        
+        @test grid_ax1_cyl == 74
+        @test grid_ax2_cyl == 1
+        @test grid_ax3_cyl == 143
+        @test grid_ax1_cart == 178
+        @test grid_ax2_cart == 178
+        @test grid_ax3_cart == 149
+    end
+    
+    # Spacing out of bounds
+    config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+    config_dict["grid"]["spacing_surface_refinement"] = [1e-2, 1e-2, 1e-2]
+    sim = @test_nowarn Simulation{T}(config_dict)
+    
+    @test_logs (
+        (:warn, r"out of bounds"s)
+    ) match_mode=:any begin
+        timed_calculate_electric_potential!(sim, verbose=false, depletion_handling=true)
+    end
+          
+    # Cylindrical grid: φ-axis not refined
+    begin
+        config_dict = SolidStateDetectors.parse_config_file(SSD_examples[:IVCIlayer])
+        
+        config_dict["grid"]["axes"]["phi"] = Dict(
+            "from" => "0",
+            "to"   => "120"
+        )
+        
+        sim_phi = @test_nowarn Simulation{T}(config_dict)
+        
+        delete!(config_dict["grid"], "spacing_surface_refinement")
+        sim_phi_noref = @test_nowarn Simulation{T}(config_dict)
+        timed_calculate_electric_potential!(sim_phi_noref, verbose = false, depletion_handling = true)
+        notref_phi_len = length(sim_phi_noref.electric_potential.grid.axes[2].ticks)
+
+        config_dict["grid"]["spacing_surface_refinement"] = [1e-3,1e-3,1e-3]
+        sim_phi_ref = @test_nowarn Simulation{T}(config_dict)
+        timed_calculate_electric_potential!(sim_phi_ref, verbose = false, depletion_handling = true)
+        ref_phi_len = length(sim_phi_ref.electric_potential.grid.axes[2].ticks)
+        
+        @test ref_phi_len == notref_phi_len
+    end
+end