For renewable wave energy to operate at grid scale, large arrays of Wave Energy Converters (WECs) need to be deployed in the ocean. Due to the hydrodynamic interactions between the individual WECs of an array, the overall power absorption and surrounding wave field will be affected, both close to the WECs (near field effects) and at large distances from their location (far field effects). Therefore, it is essential to model both the near field and far field effects of WEC arrays. It is difficult, however, to model both effects using a single numerical model that offers the desired accuracy at a reasonable computational time. The objective of this paper is to present a generic coupling methodology that will allow to model both effects accurately. The presented coupling methodology is exemplified using the mild slope wave propagation model MILDwave and the Boundary Elements Methods (BEM) solver NEMOH. NEMOH is used to model the near field effects while MILDwave is used to model the WEC array far field effects. The information between the two models is transferred using a one-way coupling. The results of the NEMOH-MILDwave coupled model are compared to the results from using only NEMOH for various test cases in uniform water depth. Additionally, the NEMOH-MILDwave coupled model is validated against available experimental wave data for a 9-WEC array. The coupling methodology proves to be a reliable numerical tool as the results demonstrate a difference between the numerical simulations results smaller than 5% and between the numerical simulations results and the experimental data ranging from 3% to 11%. The simulations are subsequently extended for a varying bathymetry, which will affect the far field effects. As a result, our coupled model proves to be a suitable numerical tool for simulating far field effects of WEC arrays for regular and irregular waves over a varying bathymetry.