We study the effect of anisotropic surface free energy on thermal grain boundary grooving in three dimensions (3D) where the transport mechanism is surface diffusion. Based on Herring's physical model (The Physics of Powder Metallurgy, McGraw-Hill, New York, 1951, p. 143), we develop a computational model in 3D and simulate the grain boundary grooving of polycrystals consisting of identical grains with a rectangular base. When the anisotropy is mild, the grooving kinetics follow power laws with exponent 1/4. The ratio of the depth of the groove pit to that of groove root is sensitive to the level of anisotropy and can deviate significantly from the isotropic case. When the anisotropy is severe, surface facets emerge and coarsen in the evolution. Severe anisotropy also induces a slight oscillation in the depth of the groove root along a grain boundary. The simulated grooving surface morphology is similar to the experimental observations of grooving with tungsten. (c) 2005 Elsevier B.V. All rights reserved.