diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index 0eb0538aa..fa594c63a 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -301,9 +301,9 @@ function add_box!(
                 Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
             end
             if segments !== nothing
-                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T, segments)
+                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T, segments; ztop = ztop, zbot = zbot)
             else
-                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
+                Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T; ztop = ztop, zbot = zbot)
             end
         end
         # Set the phase. Different routines are available for that - see below.
@@ -449,7 +449,7 @@ function add_layer!(
     if !isempty(ind_flat)
         # Compute thermal structure accordingly. See routines below for different options
         if !isnothing(T)
-            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
+            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T; ztop = maximum(zlim), zbot = minimum(zlim))
         end
 
         # Set the phase. Different routines are available for that - see below.
@@ -521,7 +521,7 @@ function add_sphere!(
     if !isempty(ind_flat)
         # Compute thermal structure accordingly. See routines below for different options
         if T != nothing
-            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
+            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T; ztop = cen[3] + radius, zbot = cen[3] - radius)
         end
 
         # Set the phase. Different routines are available for that - see below.
@@ -802,7 +802,7 @@ function add_polygon!(
     if !isempty(ind)
         # Compute thermal structure accordingly. See routines below for different options
         if T != nothing
-            Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
+            Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T; ztop = maximum(zlim_), zbot = minimum(zlim_))
         end
 
         # Set the phase. Different routines are available for that - see below.
@@ -920,9 +920,9 @@ function add_plate!(
     if !isempty(ind)
         if T != nothing
             if segments !== nothing
-                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T, segments)
+                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T, segments; ztop = maximum(zlim_), zbot = minimum(zlim_))
             else
-                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
+                Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T; ztop = maximum(zlim_), zbot = minimum(zlim_))
             end
         end
         Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
@@ -1062,7 +1062,8 @@ function add_volcano!(
     # @views Temp[ind .== false] .= 0.0
     if !isempty(ind_flat)
         if !isnothing(T)
-            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T)
+            # depth is measured downwards from the topography, so the surface is at depth = 0
+            Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T; ztop = 0.0, zbot = maximum(depth[ind]))
         end
     end
 
@@ -1254,7 +1255,7 @@ Parameters
     T = 1000
 end
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::ConstantTemp)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::ConstantTemp; kwargs...)
     Temp .= s.T
     return Temp
 end
@@ -1263,7 +1264,11 @@ end
 """
     LinearTemp(Ttop=0, Tbot=1000)
 
-Set a linear temperature structure from top to bottom
+Set a linear temperature structure from top to bottom.
+
+`Ttop` is applied at the top of the region the structure is added to (e.g. the top of the
+box or polygon) and `Tbot` at its bottom, wherever these are located; the temperature
+varies linearly in between.
 
 Parameters
 ===
@@ -1276,20 +1281,31 @@ Parameters
     Tbot = 1350
 end
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearTemp)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearTemp; ztop = nothing, zbot = nothing)
     @unpack Ttop, Tbot = s
 
-    dz = Z[end] - Z[1]
+    # Anchor at the geometric top/bottom of the region when the caller knows them;
+    # otherwise fall back to the extent of the grid points inside it.
+    zbot_grid, ztop_grid = extrema(Z)
+    ztop = isnothing(ztop) ? ztop_grid : ztop
+    zbot = isnothing(zbot) ? zbot_grid : zbot
+    dz = ztop - zbot
     dT = Tbot - Ttop
 
-    Temp = abs.((Z .- Z[end]) ./ dz) .* dT .+ Ttop
+    if dz == 0
+        return fill!(Temp, Ttop)   # region without vertical extent: no gradient to apply
+    end
+    Temp = abs.((Z .- ztop) ./ dz) .* dT .+ Ttop
     return Temp
 end
 
 """
     HalfspaceCoolingTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)
 
-Sets a halfspace temperature structure in plate
+Sets a halfspace temperature structure in plate.
+
+`Tsurface` is applied at the top of the region the structure is added to (e.g. the top of the
+box or polygon), and the depth used in the cooling model is measured from there.
 
 Parameters
 ========
@@ -1306,18 +1322,20 @@ Parameters
     Adiabat = 0        # Adiabatic gradient in K/km
 end
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp; ztop = nothing, zbot = nothing)
     @unpack Tsurface, Tmantle, Age, Adiabat = s
 
     kappa = 1.0e-6
     SecYear = 3600 * 24 * 365
-    dz = Z[end] - Z[1]
+    # Depth is measured from the geometric top of the region (the surface of the
+    # cooling halfspace), not from the topmost grid point inside it.
+    ztop = isnothing(ztop) ? maximum(Z) : ztop
     ThermalAge = Age * 1.0e6 * SecYear
 
     MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z)    # Adiabatic temperature of mantle
 
     for i in eachindex(Temp)
-        Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i] - Z[end]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i]
+        Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i] - ztop) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i]
     end
     return Temp
 end
@@ -1351,7 +1369,7 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio
     maxAge = 60       # maximum thermal age of plate [Myrs]
 end
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp; kwargs...)
     @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s
 
     kappa = 1.0e-6
@@ -1416,7 +1434,7 @@ The thermal age is capped at `maxAge` years, and the temperature is adjusted bas
 
 """
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}})
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}; kwargs...)
     @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = s
     kappa = 1.0e-6
     SecYear = 3600 * 24 * 365
@@ -1548,7 +1566,7 @@ struct Thermal_parameters{A}
     end
 end
 
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp; kwargs...)
     @unpack Tsurface, Tpot, dTadi, ubound, lbound, utbf, ltbf, age,
         dtfac, nz, rheology = s
 
@@ -1839,7 +1857,7 @@ Parameters
 - `Phase`: Phase array
 - `s`:    `McKenzie_subducting_slab`
 """
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab; kwargs...)
     @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
 
     # Thickness of the layer:
@@ -1913,7 +1931,7 @@ can be used to smooth the temperature field from continent ocean:
 - compute the thermal fields {F1} {F2}
 - then modify F.
 """
-function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)
+function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature; kwargs...)
     @unpack w_min, w_max, crit_dist, dir = s
     @unpack F1, F2 = s
 
@@ -1928,8 +1946,8 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemper
     # compute the 1D thermal structures
     Temp1 = zeros(size(Temp))
     Temp2 = zeros(size(Temp))
-    Temp1 = compute_thermal_structure(Temp1, X, Y, Z, Phase, F1)
-    Temp2 = compute_thermal_structure(Temp2, X, Y, Z, Phase, F2)
+    Temp1 = compute_thermal_structure(Temp1, X, Y, Z, Phase, F1; kwargs...)
+    Temp2 = compute_thermal_structure(Temp2, X, Y, Z, Phase, F2; kwargs...)
 
     # Compute the weights
     weight = w_min .+ (w_max - w_min) ./ (crit_dist) .* (dist)
@@ -2289,7 +2307,8 @@ function add_slab!(
 
         # Compute thermal structure accordingly. See routines below for different options {Future: introducing the length along the trench for having lateral varying properties along the trench}
         if !isnothing(T)
-            Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T)
+            # d is the distance perpendicular to the slab: its surface is at d = 0 and its base at d = -Thickness
+            Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T; ztop = 0.0, zbot = -trench.Thickness)
         end
 
         # Set the phase
diff --git a/test/test_lamem.jl b/test/test_lamem.jl
index 03f2cd623..47d83ac86 100644
--- a/test/test_lamem.jl
+++ b/test/test_lamem.jl
@@ -131,13 +131,13 @@ add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (0, 20), phase = ConstantPh
 # Linear T
 Temp = ones(Float64, size(Grid.X)) * 1350;
 add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (-50, 0), phase = ConstantPhase(3), DipAngle = 10, T = LinearTemp(Tbot = 1350, Ttop = 200))
-@test sum(Temp) == 1.1879553307069368e9
+@test sum(Temp) ≈ 1.1880350174476972e9
 
 # Halfspace cooling T structure
 Phases = zeros(Int32, size(Grid.X));
 Temp = ones(Float64, size(Grid.X)) * 1350;
 add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (-500, 0), phase = LithosphericPhases(Layers = [15 15 250], Phases = [1 2 3 0], Tlab = 1250), DipAngle = 10, T = HalfspaceCoolingTemp(Age = 20, Adiabat = 0.3))
-@test sum(Temp) == 1.194094895654067e9
+@test sum(Temp) ≈ 1.1942982365477426e9
 
 # Mid-oceanic ridge cooling temperature structure
 Phases = zeros(Int32, size(Grid.X));
diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl
index aec3164c0..fae542cbf 100644
--- a/test/test_setup_geometry.jl
+++ b/test/test_setup_geometry.jl
@@ -226,7 +226,7 @@ add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-8
 T_slab = LinearWeightedTemperature(crit_dist = 600, F1 = TsHC, F2 = TsMK);
 Temp = ones(Float64, size(Cart)) * 1350;
 add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-80.0, 0.0), phase = ConstantPhase(5), T = T_slab);
-@test sum(Temp) ≈ 3.496951166102279e8
+@test sum(Temp) ≈ 3.499457641038468e8
 
 
 Data_Final = addfield(Cart, "Temp", Temp)
@@ -248,7 +248,7 @@ add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-8
 # add accretionary prism
 add_polygon!(Phase, Temp, Cart; xlim = (500.0, 200.0, 500.0), ylim = (100.0, 400.0), zlim = (-5.0, -5.0, -65.0), phase = ConstantPhase(8), T = LinearTemp(Ttop = 20, Tbot = 30))
 @test maximum(Phase) == 8
-@test minimum(Temp) == 20.0
+@test minimum(Temp) ≈ 20.22486772486773
 @test sum(Phase) == 293264
 
 # Test the Bending slab geometry
@@ -275,7 +275,7 @@ TsHC = HalfspaceCoolingTemp(Tsurface = 20.0, Tmantle = 1350, Age = 30, Adiabat =
 temp = TsHC;
 
 add_slab!(Phase, Temp, Cart, t1, phase = phase, T = TsHC)
-@test Temp[84, 84, 110] ≈ 1042.7807110443487
+@test Temp[84, 84, 110] ≈ 1045.1322688510577
 @test extrema(Phase) == (1, 4)
 
 # with weak zone
@@ -303,7 +303,7 @@ phase = LithosphericPhases(Layers = [5 7 88], Phases = [2 3 4], Tlab = nothing)
 t1 = Trench(Start = (400.0, 400.0), End = (800.0, 800.0), θ_max = 90.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0, d_decoupling = 100.0, type_bending = :Ribe, WeakzoneThickness = 10, WeakzonePhase = 9)
 
 add_slab!(Phase, Temp, Cart, t1, phase = phase, T = T_slab)
-@test Temp[84, 84, 110] ≈ 623.9868388771819
+@test Temp[84, 84, 110] ≈ 624.6682008876219
 
 Data_Final = CartData(X, Y, Z, (Phase = Phase, Temp = Temp))
 
@@ -322,7 +322,7 @@ add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T = HalfspaceC
 T_slab = LinearWeightedTemperature(F1 = HalfspaceCoolingTemp(Age = 40), F2 = McKenzie_subducting_slab(Tsurface = 0, v_cm_yr = 4, Adiabat = 0.0), crit_dist = 600)
 add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T = T_slab);
 
-@test sum(Temp) ≈ 8.571247449404927e7
+@test sum(Temp) ≈ 8.571402268095453e7
 @test extrema(Phases) == (0, 2)
 
 # Add them to the `CartData` dataset:
@@ -365,7 +365,7 @@ add_slab!(Phases, Temp, Grid2D, trench, phase = lith, T = T_slab);
 ind = findall(Temp .> 1250 .&& (Phases .== 2 .|| Phases .== 5));
 Phases[ind] .= 0;
 
-@test sum(Temp) ≈ 8.291641108619794e7
+@test sum(Temp) ≈ 8.292000736425713e7
 @test extrema(Phases) == (0, 6)
 #Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp))
 #write_paraview(Grid2D,"Grid2D_SubductionCurvedOverriding");
@@ -465,7 +465,41 @@ add_ellipsoid!(PhasesV, TempV, Grid, cen = (4, 15, -17), axes = (1, 2, 3), Strik
 
 # Add data to cell fields:
 add_box!(PhasesC, TempC, Grid, xlim = (2, 4), zlim = (-15, -10), phase = ConstantPhase(3), DipAngle = 10, T = LinearTemp(Tbot = 1350, Ttop = 200), cell = true)
-@test sum(TempC[1, 1, :]) ≈ 13235.793377972634
+@test sum(TempC[1, 1, :]) ≈ 13051.985346405856
 
 add_ellipsoid!(PhasesC, TempC, Grid, cen = (4, 15, -17), axes = (1, 2, 3), StrikeAngle = 90, DipAngle = 45, phase = ConstantPhase(2), T = ConstantTemp(1600), cell = true)
-@test all(extrema(TempC) .≈ (200, 1600.0))
+@test all(extrema(TempC) .≈ (253.34427030131295, 1600.0))
+
+
+# ---------------------------------------------------------------------------
+# Thermal structures must be anchored at the *geometric* top of the region, not
+# at the topmost grid point inside it (which depends on the grid resolution).
+# See PR #206 (Ttop/Tsurface wrongly assumed at z = 0) and its follow-up fix.
+erfc = GeophysicalModelGenerator.erfc
+halfspace_T(depth_km, Tsurf, Tmantle, Age_Myr) = (Tsurf - Tmantle) * erfc(depth_km * 1.0e3 / (2 * sqrt(1.0e-6 * Age_Myr * 1.0e6 * 3600 * 24 * 365))) + Tmantle
+
+# 1) Box with its top at z = 0, on a grid WITH a node at z = 0 and on one WITHOUT:
+#    the temperature at the same physical depth must be the analytic halfspace value in both cases
+x = -500.0:100.0:500.0; y = -10.0:20.0:10.0
+for z in (-1000.0:20.0:50.0, range(-1000.0, 50.0, length = 65))     # 2nd grid: 16.4 km spacing, no node at 0
+    Cart = CartData(xyz_grid(x, y, z))
+    Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x))
+    add_box!(Phase, Temp, Cart; xlim = (-500.0, 500.0), zlim = (-1000.0, 0.0), phase = ConstantPhase(3), T = HalfspaceCoolingTemp(Tsurface = 0, Tmantle = 1350, Age = 100))
+    k = findlast(z .<= 0.0)                                             # topmost node inside the box
+    @test Temp[6, 1, k] ≈ halfspace_T(-z[k], 0, 1350, 100)
+    @test Temp[6, 1, k - 5] ≈ halfspace_T(-z[k - 5], 0, 1350, 100)
+end
+
+# 2) Box whose top is at depth (the case PR #206 addressed): Tsurface applies at the box top, not at z = 0
+z = range(-1000.0, 50.0, length = 65)
+Cart = CartData(xyz_grid(x, y, z))
+Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x))
+add_box!(Phase, Temp, Cart; xlim = (-500.0, 500.0), zlim = (-1000.0, -10.0), phase = ConstantPhase(3), T = HalfspaceCoolingTemp(Tsurface = 0, Tmantle = 1350, Age = 100))
+k = findlast(z .<= -10.0)
+@test Temp[6, 1, k] ≈ halfspace_T(-10.0 - z[k], 0, 1350, 100)
+
+# 3) Polygon with its top at depth and a linear gradient: Ttop at the polygon top, Tbot at its bottom
+Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x))
+add_polygon!(Phase, Temp, Cart; xlim = (-400.0, 400.0, 400.0, -400.0), ylim = (-10.0, 10.0), zlim = (-10.0, -10.0, -210.0, -210.0), phase = ConstantPhase(2), T = LinearTemp(Ttop = 20, Tbot = 420))
+k = findlast(z .<= -10.0)
+@test Temp[6, 1, k] ≈ 20 + (-10.0 - z[k]) / 200.0 * 400        # 2 C/km below the polygon top