Strictly speaking, my question is about the 1d case but I'm really interested in the 2d case so I'm leaving it open. However, I'm giving here a (partial) answer for the 1d case.

Consider the 1d problem of minimizing
$$
\int_{-1}^1 g(x) u(x) + \phi(u').
$$
All of the problems in the question can be reformulated in this form. Denote by $p = u'$, and note that $u$ can be recovered from $p$ by integration, using the boundary conditions (either from $u(-1)=0$ or $u(1)=0$). The strong form is
$$(\phi'(p(x)))' = g(x) \text{ on } (-1,1).$$
Here, it is simplest to assume that $\phi(p)$ and $g(x)$ are even functions (which is the case of examples in the original question), in which case we can somehow integrate and also reduce the domain, to arrive at:
$$\phi'(p) = \int_0^x g(t) \, dt =: G(x) \text{ for } x \in (0,1).$$
Assuming $\phi$ is strictly convex, then $\phi'$ is invertible and $p(x) = \phi'^{-1}(G(x))$ which is an "explicit" formula for the solution. In particular, if $G(x)$ is almost everywhere finite, then $p(x)$ is in the domain where $\phi$ is finite. Since the question is focused on the case where $\phi$ does indeed take on infinite values, if the domain of $\phi$ is a bounded interval, this proves that $p$ is bounded, and then $u$ is in $W^{1,\infty}$, which is already some sort of regularity. Furthermore, if $\phi'^{-1}$ and $G$ are both differentiable, then by the chain rule, so is $p$ and then $u$ is twice differentiable, etc...

We can see from the above that the case where the domain of $\phi$ is bounded is somehow *better* than the usual case where the domain of $\phi$ is $\mathbb{R}$, I am very much wondering whether that is the case in higher dimensions. I would like to believe that it is, but I can't find this stuff in the literature.