The fact of global glacier mass loss raises concerns about the sea level rise, the management of local freshwater resources, and the potential geohazards. It is imperative to understand the current glacier dynamics and to predict their evolutions. Here, we present an open source two-dimensional thermomechanically coupled ice flow model called PoLIM (i.e. Polythermal Land Ice Model). PoLIM includes higher-order stress gradients and is capable of modeling the dynamics of glaciers with higher aspect ratios. We use an enthalpy gradient method to represent the energy balance equation. This method simplifies the numerical implementation of the energy balance equation and the calculation of the cold-temperate transition surface. We also implement a drainage model to simulate the water transport in temperate ice driven by gravity and pressure gradient. A subglacial hydrology model is also included in PoLIM. We validate PoLIM against the standard benchmark experiments. It shows that PoLIM performs well in the benchmarks, indicating that its capability of simulating the thermomechanical behaviours of glaciers and subglacial water flow. PoLIM is wished being a useful and efficient tool for the glaciology community to study the issues related to glacier dynamics.