The mixing and merging of buoyant plumes originating from multiple small cooling towers into the atmosphere is numerically investigated. The effects of different arrangements of cooling towers as well as outlet geometries on the mixing of the plumes are examined. The side by side and tandem arrangements of two sources and also two types of multi-flue cooling towers are considered. The various ways by which the counter rotating vortex pair, as the dominant mechanism, affect the flow pattern in each aforementioned configuration are investigated. For tandem arrangement, far from sources, the outlet flow of the downstream cooling tower surrounds the plumes originating from the upstream cooling tower and in the region near the cooling towers, the pollutants are mostly originated from the upstream cooling tower. Maximum pollutant concentrations at distances 10 and 40 times the diameter downstream of the leading cooling tower increase by 67% and 29% with respect to those of a single cooling tower, respectively. For the side by side arrangement, the counter rotating vortices are stretched due to the large low pressure area created downstream of the cooling towers. Mixing of the plumes with the surrounding air is reduced as a result of contraction of vortices. Maximum contaminant concentrations at distances 10 and 40 times the diameter downstream of the cooling towers increase by 29% and 41% with respect to those of a single tower, respectively. Finally, the differences between flow fields formed around diamond and square configurations of multi-flue cooling towers are extensively discussed.