A 3D geological model was created for west central Alberta including 50 surfaces from the Precambrian to the bedrock topography. Uncertainty analysis for geological surfaces is used to provide information about the reliability of the 3D geological model. The uncertainty in surface modeling is a result of the error in estimation which is defined as the difference between the predicted and the observed values. The prediction of uncertainty can be measured by geostatistical tools to define the accuracy of the surface model to honor the geological points. Each geological surface was modeled using both ArcMap and Petrel. The benefit of using the kriging algorithm in ArcMap is that the interpolation result are reasonably accurate, and uncertainty analysis is easily obtained as a spatial map showing the standard deviation of the kriging estimate. However, this method tends to produce unrealistic results in areas of geological complexity. The surfaces were also modelled in Petrel using convergent interpolator which produced a more realistic representation in areas of complex geology. Unfortunately it was difficult to assess the uncertainty for the surfaces produced using this algorithm. To solve this problem, a unique workflow for assessing prediction uncertainty was developed using a combination of ArcMap, Matlab and Petrel. To get an estimate of global uncertainty, a Matlab code was written to search for the location of each observed data point and its four surrounding predicted surface points, and use their average as an estimated value. Global estimate of uncertainty (i.e. Root-mean-square error (RMSE)) do not give any information about specific areas of high and low uncertainty. However, this can be shown by producing a map of local uncertainty. The workflow developed to produce a local uncertainty map based on the convergent interpolator results starts with generating multiple subsets of the reference data set in ArcMap and then modelled in Petrel. The standard deviation of these surfaces at each location is calculated in Matlab. The results are mapped to show the cross validation of standard deviation for multiple realizations, which represents the estimation uncertainty in the surface modelled using convergent interpolator in Petrel. Uncertainty analysis leads to a better communication between the geologist and geomodelers resulting in a better understanding of the problematic areas in surface modeling.