Let $\Omega$ be a bounded domain.
If $\int_\Omega |\nabla (u-M)^+|^2 = 0$ where $u \in H^1(\Omega)$, and $u|_{\partial\Omega} \leq M$, does it follow that $u \leq M$ a.e?
I do not think this is true since we could have $u$ to be a constant bigger than $M$ except on hte boundary, which is null.
This is true, with an appropriate interpretation of "$u\le M$ on $\partial \Omega$". When working with Sobolev functions, you don't want to just define them on $\partial \Omega$ independently of what they do in $\Omega$; the resulting object would be an uninteresting combination of two unrelated functions.
Instead, the inequality on the boundary should be understood in the sense of trace operator. Without going into technicalities, having boundary values $\le M$ amounts to $(u-M)^+\in H_0^1(\Omega)$; you can take this as a definition of "$u\le M$ on $\partial \Omega$".
Then the conclusion is easy: a $H_0^1$ function with zero gradient must be identically zero. Indeed, it is constant a.e. (follows by mollification, or from the ACL property), and this constant must be $0$.