As you noted, the authors use single bars $|\cdot |$ for Hilbert space norm, without any subscript. And they use $|\cdot |$ for pointwise vector norm as well, as in (2.11). This made me look for the downvote button, but the journal's site does not have one. Anyway, the first step after the triangle inequality is Cauchy-Schwarz:
$$\|u\cdot \nabla C\|_{L^2(\Omega)}=\left(\int_\Omega |u\cdot \nabla C|^2\right)^{1/2} \le \left(\int_\Omega |u|^4\right)^{1/4}\left(\int_\Omega |\nabla C|^4\right)^{1/4}$$
The second step is the estimate
$$\|u\|_{L^4}\le M\|u\|_{L^2}^{1/2}\|\nabla u\|_{L^2}^{1/2}\tag{1}$$
(and the same for $\nabla C$). At first glance, (1) looks like the Sobolev embedding combined with trivial interpolation between $L^2$ and $L^\infty$ but it's not because $H^1$ does not embed into $L^\infty$ (in nice planar domains it embeds into $L^p$ for $p<\infty$). However, (1) is still true: this is a special case of Gagliardo-Nirenberg interpolation inequality, stated, for example, on page 314 of the book Functional Analysis, Sobolev Spaces and Partial Differential Equations by Haim Brezis. Brezis does not give a proof and refers the reader to Friedman's book which I don't have... but I just noticed it was republished by Dover! Extra special +1 to your question for this.
The Sobolev inequalities allow one to trade regularity (in the sense of derivatives) for integrability. Namely, you can trade a derivative at a lower integrability exponent $p$ for a higher integrability exponent $q=p^*=\frac{np}{n-p}>p$ (exponents $q$ in between can be obtained from Hölder's inequality if $U$ is assumed to be bounded as you do). You can iterate this inequality as many derivatives as you want, of course. Moreover, this inequality is sharp (see the last paragraph below). By looking e.g. at the proof by Gagliardo and Nirenberg reproduced e.g. in Theorem 5.6.1, pp. 277-279 of the book by Lawrence C. Evans, Partial Differential Equations, 2nd. ed., AMS, 2010, one sees that this trade comes essentially from the fundamental theorem of Calculus.
(Edit: Such an inequality is aided by another (also called Sobolev's) inequality, which is valid for $p>n$: suppose that $U$ and $u$ are as you stated for such a $p$ - then $u\in L^\infty(U)$ and satisfies the estimate $$\|u\|_{L^\infty(U)}\leq C\|u\|_{W^{1,p}(U)}\ .$$ In fact, $u$ is then even continuous since smooth elements of $W^{1,p}(U)$ are dense therein. The proof of this inequality is simple - after multiplying $u$ by a radial cutoff function $\psi$ supported in the interior $\mathrm{int}\ U$ of $U$, centered around a given point $x$ of $\mathrm{int}\ U$ and with integral 1, one passes to spherical coordinates centered at $x$ and applies the fundamental theorem of Calculus to the radial coordinate. The argument is closed by applying Hölder's inequality (for details, see e.g. the Appendix from the book by Christopher D. Sogge, Lectures on Non-Linear Wave Equations, 2nd. ed., International Press, 2008). Once more, one sees that the allowed trade between regularity and integrability comes from the fundamental theorem of Calculus. Higher derivatives are obtained by applying the Gagliardo-Nirenberg-Sobolev inequality you stated to the rhs of the above as many times as needed. Conversely, the Gagliardo-Nirenberg-Sobolev inequality for $p>1$ can be obtained from a similar argument - plus a technical lemma due to Hardy, Littlewood and Sobolev regarding the boundedness of convolution with the function $g(x)=\|x\|^{1-n}$ -, which is actually due to Sobolev himself)
What is this important for? Recall that, unlike for $p=\infty$, derivatives in Sobolev spaces for $p<\infty$ are distributional or weak derivatives: namely, if $f\in W^{1,p}$, one defines the weak partial derivative $\partial_j f$ in the $j$-th variable by formally integrating by parts against a smooth function $\phi$ compactly supported in the interior $U$ (recall that the boundary of $U$ has zero $n$-dimensional Lebesgue measure and weak derivatives need only to be defined almost everywhere): $$\int_U (\partial_j f(x))\phi(x)\mathrm{d}^nx\doteq-\int_U f(x)\partial_j\phi(x)\mathrm{d}^nx\ .$$ The definition of $W^{1,p}(U)$ states that $\partial_j f\in L^p(U)$ for all $j$.
Using Sobolev's inequality, if we set $q=\infty$ (see Edit above!) we turn the remaining, "untraded" weak derivatives of $f$ into classical derivatives. This is of extreme importance in the study of linear partial differential equations (for which Sobolev spaces were originally devised, by the way), for it is easy to obtain weak solutions in Sobolev spaces (that is, the derivatives in the partial differential operator act as weak derivatives) by using certain a priori inequalities depending on the kind of operator you are studying (elliptic, hyperbolic, etc.). If weak solutions have enough regularity, one can use the Sobolev inequality to show they are solutions in the usual sense - that is, the derivatives in the partial differential operator actually act on such solutions as classical derivatives. Check my answer to this other math.SE question of yours for an example of this.
In the case where $U$ is no longer bounded, since integrability entails a certain decay towards infinity, the bounds on the possible trade between integrability and regularity allowed by Sobolev's inequality on such an $U$ may be seen physically as an uncertainty principle in the sense of quantum mechanics - in fact, the family of uncertainty principles derived from Sobolev's inequality is far more precise than Heisenberg's original statement and plays a key role in the proof of stability of quantum matter.
Ah yes, I almost forgot to mention its important geometrical interpretation: if you take $n>1$, $f$ equal to the characteristic function of $U$, $p=1$ and $q=p^*=\frac{n}{n-1}$, the Sobolev inequality becomes the isoperimetric inequality in $\mathbb{R}^n$, which states that the isoperimetric ratio $$\frac{V^{\frac{n-1}{n}}}{A}\ ,$$ where $V$ is the $n$-dimensional volume of a bounded region $U$ with finite perimeter and $A$ is the $(n-1)$-dimensional area of its boundary (here "finite perimeter" means that $A$ is well defined and finite), is less than or equal to the same ratio for a $n$-dimensional ball (= the so-called isoperimetric constant of $\mathbb{R}^n$), with equality only in this case - that is, the inequality is sharp, as claimed in the first paragraph. The isoperimetric constant of $\mathbb{R}^n$, by the way, is the optimal value of the constant $C$ for $p=1$. The argument can be extended to other Riemannian manifolds which enjoy some form of the isoperimetric inequality (possibly with a different isoperimetric constant). For a thorough discussion on these matters, check out the excellent book Some Nonlinear Problems in Riemannian Geometry by Thierry Aubin (Springer, 1998).
Best Answer
If $p \geq 2$, then an appropriate version of the Hölder inequality shows that for $U$ bounded, $L^p(U) \hookrightarrow L^2(U)$, hence $W^{1,p}\hookrightarrow W^{1,2} \hookrightarrow L^2$ using the trace inequality. So you have
$$ \|u\|_2 \leq \|u\|_2^{1-1/p}\|u\|_2^{1/p} \lesssim \|u\|_p^{1-1/p} \|u\|^{1/p}_{W^{1,2}} \lesssim \|u\|^{1-1/p}_p \|u\|^{1/p}_{W^{1,p}} $$
which does not require passing through the Gagliardo-Nirenberg inequality. (Though one can argue that the proof of the GN inequality can be recycled to prove the trace theorem, so the two aren't completely unrelated.)
The interesting case is when $p < 2$, in which your estimate follows partially from the generalized Sobolev imbedding theorems (which allows also for trace estimates). See, for example, Theorem 4.12 Case C of Robert Adams' Sobolev Spaces. In particular the numerology requires $p \leq 2 \leq (n-1)p / (n-p)$ which implies
$$ \frac{2n}{n+1} \leq p \leq 2 $$
must hold for the classical trace theorem $W^{1,p}(\Omega) \hookrightarrow L^2(\partial\Omega)$. (If you believe in fractional Sobolev spaces, then you can also "obtain" the above by combining a fractional version of Gagliardo-Niremberg from $W^{1,p}\hookrightarrow H^s$ for some $s > 1/2$ and using the fractional trace theorem $H^s(\Omega)\hookrightarrow L^2(\partial\Omega)$; note that this route naively will require the lower bound on $p$ to be a strict inequality due to the failure of $H^{1/2}(\Omega)\not\hookrightarrow L^2(\partial\Omega)$.) In particular, for dimension $n > 1$ the $L^1$ endpoint is not covered by this theorem.
This failure at $L^1$ is not a problem of the method: it is a genuine failure of the inequality. Consider the function in two dimension $u(x,y) = \left[ (x-1)^2 + y^2\right]^{1/4}$, which is just a translated version of the function $\frac{1}{\sqrt{r}}$. It is not $L^2$ integrable when restricted to the circle, since it has a Logarithmic singularity at $(1,0)$. But $\frac{1}{\sqrt{r}}$ is certainly $L^1_{loc}(\mathbb{R}^2)$, and its derivative which has a $\frac{1}{r^{3/2}}$ singularity is also $L^1_{loc}(\mathbb{R}^2)$. (The same example also works in arbitrary $n \geq 2$; you can either use the same principle via functions like translates of $\frac{1}{r^{(n-1)/2}}$, or just argue through the method of descent.)