The shape of slabs in the mantle combined with mechanical analysis (analytic or numerical modeling) can provide information on the state of stress in the slab (e.g., compressional versus tensional), the history of subduction (e.g., rate and direction), interaction of the slab with mantle structure (e.g., layers with different material properties) or with large-scale mantle deformation driven by heating from within. Because the primary source of information about slab location is seismicity, the shapes of slabs are often reconstructed from inhomogeneously distributed and noisy observations. Previous attempts to reconstruct three-dimensional slab shapes have relied on approximating the entire slab structure using the analytic solution for a thin elastic or viscous sheet, a surface with Gaussian curvature, or smooth splines. However, these approaches enforce a long wavelength smooth shape on the slab structure that can miss shorter length-scale features, or can introduce short wavelength undulations that are unrelated to the actual slab shape. As the slab structure is often used to create input for analytic or numerical models (e.g., thermal and mechanical structure for deformation models), errors in the slab structure can lead to errors in the predicted deformation. We present a new method for generating 3D slab morphology from seismicity that generates neighborhood surface patches from algebraic surfaces of specified degree (quadratic to quintic) and then blends adjacent patches using a moving least-squares algorithm to create a smooth, continuous surface. The new procedure provides options for pre-processing noisy data to remove outliers, to systematically combine different data sets with user defined weighting, and for defining the 3D edge of the slab surface. Generation of the 3D slab surface using algebraic surface patches and blending allows for smooth, continuous assignment of properties (e.g., temperature or viscosity) that are spatially-related to the location of the slab surface and its edges (e.g., by Euclidean distance) to a numerical model grid without aliasing.