Modeling the dynamics of Greenland and Antarctic ice sheets is critical for accurately modeling climate and in particular for providing projections of sea level rise. At the core of ice-sheet modeling is the Stokes system of equations for computing the ice velocity and pressure. This is typically the most expensive part of an ice-sheet model and several less expensive lower-fidelity models have been proposed over the years to approximate the Stokes equations and enable large-scale simulations. These models rely on the fact that ice sheets are very shallow. A recently published approximation, the so-called Mono Layer High-Order model, combines a classic finite element discretization in the horizontal direction with a prescribed velocity profile (modal basis) in the (shallow) vertical direction. In this work we implement a similar model, using a tensor-product finite element on prisms, where we combine a linear nodal finite element for the triangular base of the prism with a high-order polynomial profile (the modal basis) in the vertical direction. We then select the optimal modal basis by minimizing the difference of the velocity computed with the proposed discretization with the velocity obtained with a classic, converged, finite element approximation. For accelerating the assembly we use sum-factorization techniques that expose the tensor-product structure of the finite element basis and of the quadrature rules. The effectiveness of the reduced-order model will be shown targeting the Greenland and Antarctic ice sheets.