The author's assertion
$$\underset{|p| \le C}{\sup} \int_\Omega \lambda \vec{\nabla}\cdot \vec{p}\, dx = \int_\Omega C|\nabla\lambda|\, dx \tag1$$
is true as long as $\lambda$ is sufficiently regular; being in the Sobolev class $W^{1,1}(\Omega)$ is enough. (Being Lipschitz, or having bounded continuous gradient, is more than enough.)
Inequality in one direction follows from integration by parts:
$$ \int_\Omega \lambda \vec{\nabla}\cdot \vec{p}\, dx =
- \int_\Omega \nabla \lambda \cdot \vec{p}\, dx \tag2$$ and slapping $|\cdot|$ on the right. In the opposite direction, we must produce $\vec p$ with $|\vec p|\le C$ for which $-\nabla\lambda \cdot \vec{p}$ is equal or very close to $C|\nabla\lambda|$. Formally, $\vec p=-C\nabla \lambda/|\nabla \lambda|$ should work, but some approximation may be necessary to make $\vec p$ smooth. Details can be found in the book by Guisti cited as [20] in the paper (Minimal surfaces and functions of bounded variation). There's also Wikipedia article on Bounded variation.
The supremum on the left of (1) is the standard definition of total variation of $\lambda$; it generalizes the one-dimensional concept of bounded variation. See Two definitions of "Bounded Variation Function"
The problem
$$\underset{p}{\sup} \int_\Omega \{ \lambda(\vec{\nabla}\cdot \vec{p}) - \alpha(|p| - C)\}\, dx \tag3$$
(with appropriate $\alpha$) is a relaxation of the problem in (1). Typically, $\alpha=0$ when the argument is nonpositive, and grows rapidly for positive argument. This means that you are not prevented from taking $|p|>C$, but you'll pay a fine if you do. Relaxation is a common way to approximate constrained optimization problems with unconstrained ones. However, the relaxed problem is not "equivalent" to the original one; it is merely an approximation, which may be good or bad depending on circumstances.