The Variational Boussinesq Model (VBM) for waves (Klopman et al. 2010) is based on the Hamiltonian structure of gravity surface waves. In its approximation, the fluid potential in the kinetic energy is approximated by the sum of its value at the free surface and a linear combination of vertical profiles with horizontal spatially dependent functions as coefficients. The vertical profiles are chosen a priori and determine completely the dispersive property of the model. For coastal applications, the 1D version of the model has been implemented in Finite Element with piecewise linear basis functions and has been compared with experiments from MARIN hydrodynamic laboratory for focusing wave group running above a flat bottom (Lakturov et al. 2011) and for irregular waves running above a sloping bottom (Adytia and Groesen 2011). The 2D version of the model has been derived and implemented using a pseudo-spectral method with a rectangular grid in Klopman et al. 2007, 2010. A limitation of the later implementation is a lack of flexibility when the model deals with a complicated domain such as a harbour. Here, we will present an implementation of the model in 2D Finite Element which consistent with the derivation of the model via the variational formulation. To illustrate the accuracy of wave refraction and diffraction over a complex bathymetry, the experiment of Berkhoff et al. 1982 is used to compare the FE results with measurements.