To better handle landscape heterogeneities in distributed hydrological modeling, an earlier work proposed a discretisation based on nested levels, which leads to fully unstructured modeling meshes. Upon such a discretisation, traditional numerical solutions must be adapted, especially to describe lateral flow between the unstructured mesh elements. In this paper, we illustrated the feasibility of the numeric solution of the diffusion equation, representing groundwater flow, using unstructured meshes. Thus, a two-dimensional groundwater model (BOUSS2D), adapted to convex unstructured and irregular meshes was developed. It is based on the approximation of the two-dimensional Boussinesq equation using numeric techniques suitable for non-orthogonal grids. The handling of vertical and horizontal aquifer heterogeneities is also addressed. The fluxes through the interfaces among joined mesh elements are estimated by the finite volume method and the gradient approximation method. Comparisons between the BOUSS2D predictions and analytical solutions or predictions from existing codes suggest the acceptable performance of the BOUSS2D model. These results therefore encourage the further development of hydrological models using unstructured meshes that are capable of better representing the landscape heterogeneities.