diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 62ecfe455752ed4fa329443ac9c4ac81d89c3775..03f35a1ff3727c5eef1d4044687fe5f26d902498 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -622,7 +622,6 @@ function estimate_res!(st::L1L2TVState)
             st.T(adjoint, tdata, p1) +
             (isSnabla ? zero(p1) : st.beta * u)
     facetf(hfacet, n; g, u, nablau, p1, p2, tdata) =
-        #display(n' * p2)
         # norm2 is calculated later after both facet contributions are
         # consolidated
         n' * ((isSnabla ? st.beta * st.S(zero(u), nablau) : zero(p2)) + p2)
@@ -633,7 +632,8 @@ function estimate_res!(st::L1L2TVState)
     fs = (;st.g, st.u, nablau = nabla(st.u), st.p1, st.p2, st.tdata)
     space = st.est.space
 
-    facetV = Dict{Tuple{Int, Int}, MVector{ndims_u_codomain(st), Float64}}()
+    facetv = Dict{Tuple{Int, Int}, MVector{ndims_u_codomain(st), Float64}}()
+    st.est.data .= 0.
 
     for cell in cells(mesh)
         A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
@@ -670,29 +670,33 @@ function estimate_res!(st::L1L2TVState)
             v = A[:, j] - A[:, i]
             hfacet = norm(v)
             normal = -normalize(cross(v))
-            intel = hfacet
 
-            # we don't need to recompute values since only piecewise
+            # we don't need to recompute `opvalues` since only piecewise
             # constant data is used in `facetf`
-            # we use sqrt here since this value will be squared later
-            facetres = (isSnabla ? sqrt(hfacet) : inv(sqrt(hfacet))) *
-                facetf(hfacet, normal; opvalues...) .* sqrt(intel)
+            facetres = facetf(hfacet, normal; opvalues...)
 
             # average facet contributions
-            v = SVector(facetres)
-            facetV[(fi, fj)] = get(facetV, (fi, fj), zero(v)) + SVector(v)
-            facetV[(fj, fi)] = get(facetV, (fj, fi), zero(v)) + SVector(v)
+            v = SVector{ndims_u_codomain(st)}(facetres / 2)
+            facetv[(fi, fj)] = get(facetv, (fi, fj), zero(v)) + v
+            facetv[(fj, fi)] = get(facetv, (fj, fi), zero(v)) + v
         end
     end
 
     # add facet contributions to cells
     for cell in cells(mesh)
         for (i, j) in ((i, mod1(i + 1, 3)) for i in 1:3)
+            A = SArray{Tuple{ndims_space(mesh), nvertices_cell(mesh)}}(
+                view(mesh.vertices, :, view(mesh.cells, :, cell)))
+            v = A[:, j] - A[:, i]
+            hfacet = norm(v)
+            intel = hfacet
+
             fi = mesh.cells[i, cell]
             fj = mesh.cells[j, cell]
 
             gdofs = space.dofmap[:, 1, cell]
-            st.est.data[gdofs] .+= norm2(facetV[(fi, fj)])
+            st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) *
+                norm2(facetv[(fi, fj)]) * intel
         end
     end