Summary: A three dimensional numerical analysis for unconfined seepage problems in inhomogeneous and anisotropic domains with arbitrary geometry is presented in this paper. The unconfined seepage problems are nonlinear in its nature due to unknown location of the phreatic surface and nonlinear boundary conditions which complicates its solution. The presented method is based on the application of non-boundary-fitted meshes and is an extension of the recently proposed two dimensional smoothed fixed grid finite element method. The main objective of using this method is to facilitate solution of variable domain problems and improve the accuracy of the formulation of the boundary intersecting elements. In this method, the gradient smoothing technique is used to obtain the element matrices. This technique simplifies the solution significantly by reducing the volume integrals over the elements into area integrals on the faces of smoothing cells. To locate the free surface, an initial guess for the unknown geometry is selected and modified in each iteration to eventually satisfy nonlinear boundary condition. The application of the proposed technique for three dimensional seepage problems is carried out for different examples including rectangular, trapezoidal and semi-cylindrical dams and the results are compared with those available in the literature.