High precision geoid determination is a challenging task at the national scale. Many efforts have been conducted to determine precise geoid, locally or globally. Geoid models have different precision depending on the type of information and the strategy employed when calculating the models. This contribution addresses the challenging problem of combining different regional and global geoid models, possibly combined with the geometric geoid derived from GNSS/leveling observations. The ultimate goal of this combination is to improve the precision of the combined model. We employ fitting an appropriate geometric surface to the geoid heights and estimating its (co)variance components. The proposed functional model uses the least squares 2D bi-cubic spline approximation (LS-BICSA) theory, which approximates the geoid model using a 2D spline surface fitted to an arbitrary set of data points in the region. The spline surface consists of third- order polynomial pieces that are smoothly connected together, imposing some continuity conditions at their boundaries. In addition, the least-squares variance component estimation (LS- VCE) is used to estimate precise weights and correlation among different models. We apply this strategy to the combined adjustment of the high-degree global gravitational model EIGEN-6C4, the regional geoid model IRG2016, and the Iranian geometric geoid derived from GNSS/leveling data. The accuracy of the constructed surface is investigated with five randomly selected subsamples of check points. The optimal combination of the two geoid models along with the GNSS/leveling data shows a reduction of 21 mm (~20%) in the RMSE values of discrepancies at the check points.