Knowledge of the wave perturbation caused by an array of Wave Energy Converters (WEC) is of great concern, in particular for estimating the interaction effects between the various WECs and determining the modification of the wave field at the scale of the array, as well as possible influence on the hydrodynamic conditions in the surroundings. A better knowledge of these interactions will also allow a more efficient layout for future WEC farms. The present work focuses on the interactions of waves with several WECs in an array. Within linear wave theory and in frequency domain, we propose a methodology based on the use of a BEM (Boundary Element Method) model (namely Aquaplus) to solve the radiation-diffraction problem locally around each WEC, and to combine it with a model based on the mild slope equation at the scale of the array. The latter model (ARTEMIS software) solves the Berkhoff’s equation in 2DH domains (2 dimensional code with a z-dependence), considering irregular bathymetries. In fact, the Kochin function (a far field approximation) is used to propagate the perturbations computed by Aquaplus into Artemis, which is well adapted for a circular wave representing the perturbation of an oscillating body. This approximation implies that the method is only suitable for well separated devices. A main advantage of this coupling technique is that Artemis can deal with variable bathymetry. It is important when the wave farm is in shallow water or in nearshore areas. The methodology used for coupling the two models, with the underlying assumptions is detailed first. Validations test-cases are then carried out with simple bodies (namely heaving vertical cylinders) to assess the accuracy and efficiency of the coupling scheme. These tests also allow to analyze and to quantify the magnitude of the interactions between the WECs inside the array.