Ah, I was misled by the fact that $E\subseteq \mathbb{R}^2$. If $E$ has sufficiently smooth boundary, the usual extension and density theorems allows us to consider instead $C^\infty_0(\mathbb{R}^2)$ functions instead. Then the inequality is a simple consequence of integration by parts:
Let $f\in C^\infty_0(\mathbb{R}^d)$, then first we have the interpolation inequality
$ \int |f|^4 \mathrm{d}x \leq \left(\sup |f|\right)^2 \cdot \int |f|^2 \mathrm{d}x $
which implies
$ \|f\|_4^2 \leq \|f\|_\infty \|f\|_2 \tag{*}$
Next for the first derivative terms, let $\partial_i = \partial_{x^i}$ for some fixed coordinate index $i$, we have the following integration by parts formula
$ \int |\partial_i f|^4 \mathrm{d}x = \int |\partial_i f|^2 \partial_i f \partial_i f \mathrm{d}x = - \int f \partial_i\left( |\partial_i f|^2 \partial_i f\right) \mathrm{d}x $
So
$ \int |\partial_i f|^4 \mathrm{d}x \leq 3 \|f\|_\infty \int |\partial_i f|^2 \cdot |\partial^2_{ii} f| \mathrm{d}x \leq 3 \|f\|_\infty \| (\partial_i f)^2\|_2 \| \partial^2_{ii} f\|_2 \tag{#} $
where we used Cauchy-Schwarz in the last step. Noting that
$ \| (\partial_i f)^2\|_2 = \| \partial_i f\|_4^2 $
and that the left hand side of (#) is $\|\partial_i f\|_4^4$ we have that
$ \| \partial_i f \|_4^2 \leq 3 \|f\|_\infty \|\partial^2_{ii} f\|_2 \tag{%}$
Now summing (%) over the index $i$ and combining with (*) gives the desired expression. Notice that the inequality holds in all dimensions, not just in 2 dimensions. (The scaling behaviour of the inequality is compatible with it being dimension independent.)