Summary. In this paper, new natural element approximations are proposed, in order to address issues associated with incompressibility as well as to increase the accuracy in the Natural Element Method (NEM). The NEM exhibits attractive features such as interpolant shape functions or auto-adaptive domain of influence, which alleviates some of the most common difficulties in meshless methods. Nevertheless, the shape functions can only reproduce linear polynomials, and in contrast to moving least squares methods, it is not easy to define interpolations with arbitrary approximation consistency. In order to treat mechanical models involving incompressible media in the framework of mixed formulations, the associated functional approximations must satisfy the well known inf-sup, or LBB condition. The first proposed approach constructs richer NEM approximation schemes by means of bubbles associated with the topological entities of the underlying Delaunay tessellation, allowing to pass the LBB and to remove pressure oscillations in the incompressible limit. Despite of its simplicity, this approach does not construct approximation with higher order consistency. The second part of the paper deals with a discussion on the construction of second-order accurate NEM approximations. For this purpose, two techniques are investigated : (a) the enrichment in the MLS framework of the bubbles with higher-order polynomials and (b) the use of a new Hermite-NEM formulation.