The accurate modeling and simulation of the spread of contaminants within water distribution networks (WDNs) is an important task for drinking water security. In commonly used water quality simulation platforms the mixing of concentrations at junctions is assumed to be a complete mixing behavior. Experimental investigations have shown that this assumption is only true for certain flow conditions. Indeed, mixing depends on the geometry and the load configuration and often it is incomplete or there is some shortcut and preferential flow pathway with no mixing at all. In this paper, we present a more representative model of mixing at nodes that is valid for double T- and cross-junctions. Statistics on real WDNs are used to define realistic junction scenarios. From the latter, two-dimensional and three-dimensional computational fluid dynamics (CFD) simulations for the mixing process at different types of junctions were made and compared to experiments. Both the simulations and experiments were in agreement and show a difference in mixing of 10% and more compared to perfect mixing models. CFD results were used to build a lookup table and Kriging interpolation was applied for entries not in the table.