Gradient of Sobolev function on level sets
I know a proof which proceeds by proving Stampacchia's lemma in two different ways. I do not know if there is a more direct proof.
Therefore, we need a smoothing of the $\max(\cdot,0)$-function: $$\varphi_\varepsilon(t) = \begin{cases}0 & \text{if } t \le 0, \\ \frac{t^2}{2\,\varepsilon} & \text{if } 0 < t \le \varepsilon, \\ t - \frac\varepsilon2 & \text{if } \varepsilon < t.\end{cases}$$ (I hope I got the constants right.)
Then, you consider $g_\varepsilon := \varphi_\varepsilon \circ f$ and show that $g_\varepsilon \to \max(f,0)$ in $W^{1,p}(\Omega)$. This shows $$\nabla (\max(f,0))(x) = \begin{cases} 0 & \text{if } f(x) \le 0, \\ \nabla f(x) & \text{if } f(x) > 0.\end{cases}$$ Here, it is essential that $\varphi_\varepsilon'(0) = 0$.
Now, you perform the same calculation with the translates $\psi_\varepsilon(t) := \varphi_\varepsilon(t + \varepsilon)$. Since $\psi_\varepsilon'(0) = 1$, you end up with $$\nabla (\max(f,0))(x) = \begin{cases} 0 & \text{if } f(x) < 0, \\ \nabla f(x) & \text{if } f(x) \ge 0.\end{cases}$$
This shows that $\nabla f = 0$ a.e. on the set $\{f = 0\}$.