For 3D surface reconstruction problems with noisy and incomplete range data measured from complex scenes with arbitrary topologies, a low-level representation, such as level set surfaces, is used. Such surface reconstruction is typically accomplished by minimizing a weighted sum of data-model discrepancy and model smoothness terms. This paper introduces a new nonlinear model smoothness term for surface reconstruction based on variations of the surface normals. A direct solution requires solving a fourth-order partial differential equation (PDE), which is very difficult with conventional numerical techniques. Our solution is based on processing the normals separately from the surface, which allows us to separate the problem into two second-order PDEs. The proposed method can smooth complex, noisy surfaces, while preserving sharp, geometric features, and it is a natural generalization of edge-preserving methods in image processing, such as anisotropic diffusion.