In this article, a solution to nonlinear moving boundary problems in heterogeneous geological media using the meshless method is proposed. The free surface flow is a moving boundary problem governed by Laplace equation but has nonlinear boundary conditions. We adopt the collocation Trefftz method (CTM) to approximate the solution using Trefftz base functions, satisfying the Laplace equation. An iterative scheme in conjunction with the CTM for finding the phreatic line with over–specified nonlinear moving boundary conditions is developed. To deal with flow in the layered heterogeneous soil, the domain decomposition method is used so that the hydraulic conductivity in each subdomain can be different. The method proposed in this study is verified by several numerical examples. The results indicate the advantages of the collocation meshless method such as high accuracy and that only the surface of the problem domain needs to be discretized. Moreover, it is advantageous for solving nonlinear moving boundary problems with heterogeneity with extreme contrasts in the permeability coefficient.