Ice sheet modeling is an important tool both for reconstructing past ice sheets and predicting their future evolution, but is complex and computationally costly. It involves modeling a system including the ice sheet, ice shelves and ice streams, which all have different dynamical behavior. The governing equations are non-linear, and to capture a full glacial cycle more than 100,000 years need to be simulated. To reduce the problem size, approximations of the equations are introduced. The most common approximation, the Shallow Ice Approximation (SIA), works well in the ice bulk but fails in e.g. the modeling of ice streams and the ice sheet/ice shelf coupling. In recent years more accurate models, so-called higher order models, have been constructed to address these problems. However, these models are generally constructed in an ad hoc fashion, lacking rigor. In this thesis, so-called Second Order Shallow Ice Approximation (SOSIA) equations for pressure, vertical shear stress and velocity are implemented into the ice sheet model SICOPOLIS. The SOSIA is a rigorous model derived by Baral in 1999 [3]. The numerical solution for a simple model problem is compared to an analytical solution, and benchmark experiments, comparing the model to other higher order models, are carried out. The numerical and analytical solution agree well, but the results regarding vertical shear stress and velocity differ from other models. It is concluded that there are problems with the model implemented, most likely in the treatment of the relation between stress and strain rate.