Ice flow models used to project the mass balance of ice sheets in Greenland and Antarctica usually rely on the Shallow Ice Approximation (SIA) and the Shallow‐Shelf Approximation (SSA), sometimes combined into so‐called “hybrid” models. Such models, while computationally efficient, are based on a simplified set of physical assumptions about the mechanical regime of the ice flow, which does not uniformly apply everywhere on the ice sheet/ice shelf system, especially near grounding lines, where rapid changes are taking place at present. Here, we present a new thermomechanical finite element model of ice flow named ISSM (Ice Sheet System Model) that includes higher‐order stresses, high spatial resolution capability and data assimilation techniques to better capture ice dynamics and produce realistic simulations of ice sheet flow at the continental scale. ISSM includes several approximations of the momentum balance equations, ranging from the two‐dimensional SSA to the three‐dimensional full‐Stokes formulation. It also relies on a massively parallelized architecture and state‐of‐the‐art scalable tools. ISSM employs data assimilation techniques, at all levels of approximation of the momentum balance equations, to infer basal drag at the ice‐bed interface from satellite radar interferometry‐derived observations of ice motion. Following a validation of ISSM with standard benchmarks, we present a demonstration of its capability in the case of the Greenland Ice Sheet. We show ISSM is able to simulate the ice flow of an entire ice sheet realistically at a high spatial resolution, with higher‐order physics, thereby providing a pathway for improving projections of ice sheet evolution in a warming climate.