How to improve accuracy of enforced integral constraints involving derivatives

This is expected behavior though. With refinement, the integral of ∫∇u⋅n ds limits to h. The fact that this is not an identity relates to the operation ∫∇u⋅n ds is mathematically unbound in the relevant function space H^1, and one should not interpret that integral strongly. Instead, the definition of of the operation ∇\phi⋅n on the interface for any \phi \in H^1 and \Delta \phi \in H^{-1} actually is through intergration by parts:

∫∇u⋅n v ds := ∫∇u⋅∇v d\Omega + ∫\Delta u v d \Omega , which follows precisely from your weak form as = ∫∇u⋅∇v d\Omega + ∫f v d \Omega

There is a nice chapter “Weak normal derivatives, normal and tangential traces, and tangential differential operators on Lipschitz boundaries” by Sayas from 2009 (for which I can’t seem to find a digital link a.t.m.) explaining this. At the bottom of this pdf it is also addressed: https://www.ltcc.ac.uk/media/london-taught-course-centre/documents/Applied-Computational-Methods---Notes-4.pdf

See also this posts Consistency between displacement-controlled and load-controlled FEAs - #2 by Stein

Summarizing, in FEM, you really should not be trying to set ∫∇u^h⋅n in a strong way (it is ill-defined!). You can be perfectly happy with setting the requirement on ∫∇u⋅∇v d\Omega + ∫f v d \Omega instead, which is its mathematically correct definition.