After looking more carefully into this, I think that the fact that \nabla^2 is not defined at the boundary is not the cause of this issue.
In fact, if I enforce the BC as @Stein suggested by replacing - n[i] * (v.dx( i )) * nu_w * ds with - n[i] * gradv_expression * nu_w * ds and plot the ‘residual’ of the solution \nabla^2 w - f by setting a different color scale, it significantly differs from zero also in the bulk, where it should not:
On the other hand, if I enforce the Dirichlet BCs for w by setting
w_profile = Expression( '8 * (sin( 2 * x[0] ) - sin( 2 * x[1] ))', L=L, h=h, element=Q.sub( 2 ).ufl_element() )
bc_w = DirichletBC( Q.sub( 2 ), w_profile, boundary )
bcs = [bc_u, bc_v, bc_w]
the ‘residual’ looks fine
So there is something wrong here. Given that this is issue somewhat departs from the original question of this post, where I saw that the residual took huge values, and which have been fixed by this reply, I will open a new post where I clearly re-formulat the question about this latter point.