A numerical model is presented for the estimation of Wave Energy Converter (WEC) performance in variable bathymetry regions, taking into account the interaction of the floating units with the bottom topography. The proposed method is based on a coupled-mode model for the propagation of the water waves over the general bottom topography, in combination with a Boundary Element Method for the treatment of the diffraction/radiation problems and the evaluation of the flow details on the local scale of the energy absorbers. An important feature of the present method is that it is free of mild bottom slope assumptions and restrictions and it is able to resolve the 3D wave field all over the water column, in variable bathymetry regions including the interactions of floating bodies of general shape. Numerical results are presented concerning the wave field and the power output of a single device in inhomogeneous environment, focusing on the effect of the shape of the floater. Extensions of the method to treat the WEC arrays in variable bathymetry regions are also presented and discussed.modes may have a significant impact on the wave phase evolution during propagation. Such a fact was demonstrated through the interference process in one-directional wave propagation as observed for either varying topographies (see e.g., [12,13]) or abrupt bathymetries including coastal structures (see e.g., [14][15][16]). For such problems, the consistent coupled-mode theory has been developed in [17], for the water waves propagation in variable bathymetry regions. Furthermore, it was subsequently extended for 3D bathymetry in [18], and applied successfully to treat the wave transformation over nearshore/coastal sites with steep 3D bottom features, like underwater canyons; see, e.g., [19,20].In recent works [21,22] the coupled-mode model is further extended to treat the wave-current-seabed interaction problem, with application to the wave scattering by non-homogeneous current over general bottom topography. The problem of the directional spectrum transformation of an incident wave system over a region of strongly varying three-dimensional bottom topography is further studied in [22]. The accuracy and efficiency of the coupled-mode method is tested, comparing numerical predictions against experimental data by [23] and calculations by the phase-averaged model SWAN [24,25]. Results are shown in various representative test cases demonstrating the importance of the first evanescent modes and the additional sloping-bottom mode when the bottom slope is not negligible.In this work, a methodology is presented to treat the propagation-diffraction-radiation problem locally around each WEC, supporting the calculation of the interaction effects of the floating units with variable bottom topography at a local scale. The method is based on the coupled-mode model developed by [17], and extended to 3D by [18], for water wave propagation over general bottom topography, in conjunction with the Boundary Element Method (BEM) for the hydrodynamic analysis of fl...