For a given bounded domain \(\Omega \subseteq \mathbb{R}^N\) with a \(C^2\)-boundary \(\partial\Omega\), the aim of this paper is the study of the following nonlinear elliptic inclusion with Robin boundary condition \[ \begin{cases} -\Delta_p u(z)+\partial \varphi(u(z)) \ni f(z,u(z),\nabla u(z)) \quad & \text{in }\Omega,\\ \frac{\partial u}{\partial n_p}+\beta(z) |u|^{p-2}u=0 \quad & \text{on }\partial\Omega, \end{cases} \] where \(\Delta_p\) denotes the \(p\)-Laplace differential operator for \(2\leq p<\infty\), \(\varphi \in \Gamma_0(\mathbb{R})\), where \(\Gamma_0(\mathbb{R})\) denotes the cone of all functions \(\varphi: \mathbb{R}\to \overline{\mathbb{R}}=\mathbb{R}\cup \{\infty\}\) which are convex, lower semicontinuous and proper. Moreover, \(\partial \varphi(\cdot)\) denotes the subdifferential in the sense of convex analysis and the function \(f: \Omega \times \mathbb{R} \times \mathbb{R}^N \to \mathbb{R}\) is measurable in the first argument and locally Hölder in the others. Since \(f\) depends on the gradient of the solution (so-called convection term), variational methods cannot be applied. Therefore, by applying a topological approach based on fixed point theory (the Leray-Schauder alternative principle) and approximating the original problem using the Moreau-Yosida approximations of the subdifferential term, the authors prove the existence of a smooth solution.