The aim of this paper is to develop a new framework of surface functional models for surface functional data which contains repeated observations in both the domains. The primary interest in our problem is to investigate the relationship between a response and the two domains, where the numbers of observations in both domains within a subject may be diverging. We estimate the mean function based on local linear smoothers. Unprecedented complexity presented in the surface functional models, such as possibly distinctive sampling designs and the dependence between the two domains, makes the theoretical investigation challenging. We are able to provide a comprehensive investigation of the asymptotic properties of the mean function estimator based on a general weighing scheme. Moreover, we can mathematically categorize the surface data into nine cases according to the sampling designs of both the domains, essentially based on the relative order of the number of observations in each domain to the sample size. We derive the specific asymptotic theories and optimal bandwidth orders in each of the nine sampling design cases under all the three weighing schemes. We also examine the finite-sample performance of the estimators through simulation studies.