Continental ice sheet simulations are commonly based on either the Full Stokes (FS) model, or its simplification, the Shallow Ice Approximation (SIA) model. This thesis examines a machine learning error estimation approach for assessing the accuracy of the solutions to the SIA model, where the reference (exact) solution is that of the Stokes model. We use Gaussian Process (GP) regression through existing GP libraries in Python to model and train a machine learning model. For computational efficiency reasons we use Variational Nearest Neighbor Gaussian Processes (VNNGP), where the input data are the SIA solution and the ice sheet geometry characteristics. The output data is the error between the SIA solution and the FS solution. We find that these models trained on various ice sheet geometries are able to make rough predictions for other simple geometries not trained for; however we observe a poor fit for the much more complex Greenland geometry, which suggests further work to be done, utilizing more diverse geometries for training.