Interfaces of rather different natures---as e.g.\ bacterial-colony or forest-fire boundaries, or semiconductor layers grown by different methods (MBE, sputtering, etc)---are self-affine fractals, and feature scaling with \emph{universal} exponents (depending on the substrate's dimensionality $d$ and global topology, as well as on the driving randomness' spatial and temporal correlations but \emph{not} on the underlying mechanisms). Adding lateral growth as an essential (nonequilibrium) ingredient to the known equilibrium ones (randomness and interface relaxation), the Kardar--Parisi--Zhang (KPZ) equation succeeded in finding (via dynamic renormalization group) the correct exponents for flat $d=1$ substrates and (spatially and temporally) uncorrelated randomness. It is this \emph{interplay} which gives rise to the unique, non-Gaussian scaling properties characteristic of the specific, universal type of nonequilibrium roughening. Later on, the asymptotic statistics of process $h(x)$ fluctuations in the scaling regime was also analytically found for $d=1$ substrates. For $d>1$ substrates however, one has to rely on numerical simulations. Here we review a variational approach that allows analytical progress regardless of substrate dimensionality. After reviewing our previous numerical results in $d=1$, $2$ and $3$ on the time evolution of one of the functionals---which we call the \emph{nonequilibrium potential} (NEP)---as well as its scaling behavior with the nonlinearity parameter $\lambda$, we discuss the stochastic thermodynamics of the roughening process and the memory of process $h(x)$ in KPZ and in the related Golubovi\'c--Bruinsma (GB) model, providing numerical evidence on the significant dependence on initial conditions of the NEP's asymptotic behavior in both models. Finally, we highlight some open questions.