This paper presents a novel nonlinear hyperspectral mixture model and its associated supervised unmixing algorithm. The model assumes a linear mixing model corrupted by an additive term which accounts for multiple scattering nonlinearities (NL). The proposed model generalizes bilinear models by taking into account higher order interaction terms. The inference of the abundances and nonlinearity coefficients of this model is formulated as a convex optimization problem suitable for fast estimation algorithms. This formulation accounts for constraints such as the sum-to-one and nonnegativity of the abundances, the non-negativity of the nonlinearity coefficients, and the spatial sparseness of the residuals. The resulting convex problem is solved using the alternating direction method of multipliers (ADMM) whose convergence is ensured theoretically. The proposed mixture model and its unmixing algorithm are validated on both synthetic and real images showing competitive results regarding the quality of the inference and the computational complexity when compared to the state-of-the-art algorithms.