In this study, a meshless simulation model based on the element-free Galerkin method (EFGM) is proposed for the coupled simulation of groundwater flow and contaminant transport (GFCT) in the two-dimensions. EFGM discretizes the governing GFCT equations using the Galerkin integral approach and shape functions generated through the moving least squares (MLS) approximation method. It results into the stable system of equations for estimating the groundwater head and contaminant concentration values. EFGM computes the several global matrices using the spatial information of the scattered and unconnected field nodes and hence, terminates the meshing procedure unlike the finite difference (FDM) and finite element (FEM) methods. The Kronecker delta function property is missing in EFGM shape functions and hence, the Dirichlet boundary conditions are enforced with some special methods. We have used the penalty method to include them into the presented EFGM model. This method is straightforward and computationally useful for the real field simulations. The proposed model is first validated using both the one- and two-dimensional analytical solutions. The developed coupled model is then applied to solve the hypothetical and real field problems. EFGM results for both the problems are compared with the MODFLOW and MT3DMS results for the head and concentration values respectively. EFGM results are found to be satisfactory for both the cases, and it shows the potential of the proposed approach for the coupled GFCT simulations.