In this paper, we combine Bochner formula, Saloff-Coste’s Sobolev inequality and the Nash-Moser iteration method to study the local and global behaviors of solutions to the nonlinear elliptic equation \(\Delta _pu+\Delta _qu+h(u,|\nabla u|^2)=0\) defined on a complete Riemannian manifold \(\left( M,g\right) \) , where \(q\ge p>1\) , \(h\in C^1(\mathbb {R}\times \mathbb {R}^{+})\) and \(\Delta _z u=\textrm{div}\left( \left| \nabla u\right| ^{z-2}\nabla u\right) \) , with \(z\in \{ p,~q\}\) , is the usual z-Laplace operator. Under some assumptions on p, q and h(x, y), we derive concise gradient estimates for solutions to the above equation and then obtain some Liouville type theorems. In particular, we use integral estimate method to show that, if u is a non-negative entire solution to \(\Delta _p u +\Delta _q u=0\) ( \(n\le p\le q\) ) on a complete non-compact Riemannian manifold M with non-negative Ricci curvature and \(\dim M = n\ge 2\) , then u is a trivial constant solution.