We present a hybrid solution strategy for the numerical solution of the two-dimensional (2D) partial differential equations of Green-Nagdhi (GN), which simulates fully nonlinear, weakly dispersive free surface waves. We re-write the standard form of the equations by splitting the original system in its elliptic and hyperbolic parts, through the definition of a new variable, accounting for the dispersive effects and having the role of a non-hydrostatic pressure gradient in the shallow water equations. We consider a two-step solution procedure. In the first step we compute a source term by inverting the elliptic coercive operator associated to the dispersive effects; then in a hyperbolic step we evolve the flow variables by using the non-linear shallow water equations, with all non-hydrostatic effects accounted by the source computed in the elliptic phase. The advantages of this procedure are firstly that the GN equations are used for propagation and shoaling, while locally reverting to the nonlinear shallow water equations to model energy dissipation in breaking regions. Secondly and from the numerical point of view, this strategy allows each step to be solved with an appropriate numerical method on arbitrary unstructured meshes. We propose a hybrid finite element (FE) finite volume (FV) scheme, where the elliptic part of the system is discretized by means of the continuous Galerkin FE method and the hyperbolic part is discretized using a third-order node-centred finite volume (FV) technique. The performance of the numerical model obtained is extensively validated against experimental measurements from a series of relevant benchmark problems.