The flow around four infinite cylinders in a square arrangement was numerically modelled using dynamic Smagorinsky Large Eddy Simulation (LES). For the simulations five different pitch (P) to diameter (D) spacing ratios (P/D) were investigated at a Reynolds number of 3000. The simulated P/D ratios investigated were P/D = 1.25, 1.4, 1.5, 1.75 and 2.0. In addition to the phenomena detected for the case of two cylinders (‘‘tandem’’ and ‘‘side by side’’ in Afgan et al. (2011)), the detailed results showed the presence of a variety of complex phenomena such as jet impinging on curved surfaces with coupling layers around the cylinders. The topology of the flow changed with the change in the spacing ratios between the cylinders showing a clear biased flow behavior for spacing ratios (P/D) 1.25, 1.4 and 1.5. This bistability phenomenon was detected based on the analysis of the lift signals and the flow streamlines. Good agreement was found between the current numerical results and existing experimental data for both the local (velocity and pressure fields) and global quantities for all the simulated configurations.