We present an approach for large-scale modeling of parametric surfaces using spherical harmonics (SHs). A standard least square fitting (LSF) method for SH expansion is not scalable and cannot accurately model large 3D surfaces. We propose an iterative residual fitting (IRF) algorithm, and demonstrate its effectiveness and scalability in creating accurate SH models for large 3D surfaces. These large-scale and accurate parametric models can be used in many applications in computer vision, graphics, and biomedical imaging. As a simple extension of LSF, IRF is very easy to implement and requires few machine resources.