We present a fast method that adaptively approximates large-scale functional scattered data sets with hierarchical B-splines. The scheme is memory efficient, easy to implement and produces smooth surfaces. It combines adaptive clustering based on quadtrees with piecewise polynomial least squares approximations. The resulting surface components are locally approximated by a smooth B-spline surface obtained by knot removal. Residuals are computed with respect to this surface approximation, determining the clusters that need to be recursively refined, in order to satisfy a prescribed error bound. We provide numerical results for two terrain data sets, demonstrating that our algorithm works efficiently and accurate for large data sets with highly non-uniform sampling densities.