In this work, we present a finite element model to approximate the modified Boussinesq equations. The objective is to deal with the major problem associated with this system of equations, namely, the need to use stable velocity‐depth interpolations, which can be overcome by the use of a stabilization technique. The one described in this paper is based on the splitting of the unknowns into their finite element component and the remainder, which we call the subgrid scale. We also discuss the treatment of high‐order derivatives of the mathematical model and describe the time integration scheme.