The determination of point heights from a physical reference surface is still a challenge in geodetic and surveying applications today. Precise geoid models are used to obtain the physical point heights from the GNSS observations with height transformation. The precision of the geoid model mainly depends on the employed calculation method of the geoid as well as the accuracy, density, and distribution of the used data. This study provides a comprehensive comparison of four surface interpolation methods having different mathematical backgrounds in local geoid modeling, in the west of Turkey. In the tests, high accuracy GNSS/leveling data at the scattered control benchmarks on topography were used. The tested algorithms are multivariable polynomial regression (MPR) analysis with least-squares adjustment (LSA), stochastic based least-squares collocation (LSC), finite elements based bivariate (BIVAR) interpolation, learning-based wavelet neural networks (WNN) methods. Among these algorithms, BIVAR was applied for the first time in this study for local geoid modeling. Apart from its superiority for accuracy, this finite elements based algorithm is considerably successful than the other three methods in providing continuity of the surface model. The surface continuity is a crucial issue in local geoid modeling. In this regard, the results of this study make a significant contribution to the practical use of local geoids. The calculated local geoid models with interpolation techniques were also compared with a gravimetric geoid model in the area. In overall results, the BIVAR interpolation algorithm provided superior performance than the other tested geoid models. Hence, its use was highly recommended for local geoid modeling. The accuracy of the local geoid model calculated with the BIVAR method is 2.65 cm.