are amenable to calculation. The electrical properties can be predicted using semi-classical Boltzmann transport theory, and it has also recently become possible to use first-principles modelling to investigate the structural dynamics and lattice thermal conductivity. The latter calculations have played an important role in understanding the intrinsic anharmonicity in materials with low 𝜅 latt , including several flagship thermoelectrics, and have also been used to identify and characterize the thermal transport in novel candidate TEs. There are several reports in the literature of modelling studies on the thermoelectric properties of bulk bismuth chalcogenides and nanostructures (e.g. single layers), including a handful of studies that investigate the structural dynamics and thermal conductivity. Modelling studies on the alloys are more scarce, owing to the higher complexity of the calculations, and tend to be restricted to compositions with well-defined structures. Previous theoretical studies of alloy systems using first-principles density functional theory (DFT) and statistical thermodynamics have demonstrated the capacity systematically to explore a range of compositions between the endpoints of a two-component alloy system. The structure, properties and energetic and dynamical stability of different compositions can be evaluated in silico, providing a means to predict the electronic structure and thermal transport of the alloy and to guide experimental design towards optimizing the properties for thermoelectric applications. The high computational resource requirements for modelling the lattice dynamics and thermal transport, however, pose significant challenges to studying alloy systems this way, which has resulted in few investigations in this direction. In this work, we use first-principles DFT calculations to perform a comprehensive study of the lattice dynamics, energetic and dynamical stability and thermal transport of all of the reported phases of Bi2S3, Bi2Se3 and Bi2Te3, and we employ methodology developed in previous studies to investigate the thermodynamics and properties of Bi2(Se,S)3 and Bi2(Se,Te)3 alloys. The differing activities of the Bi cation lone pairs in Bi2S3, Bi2Se3 and Bi2Te3 lead to competition between R3 ̅ m and Pnma phases with very different thermal transport properties, which we attribute quantitatively to differences in the phonon group velocities and lifetimes. Whereas the Bi2(Se,S)3 alloy in the Pnma phase forms a homogenous solid solution, the Bi2(Se,S)3 and Bi2(Se,Te)3 alloys in the R3 ̅ m phase show strong structural preferences for certain chalcogen arrangements. This leads to significant deviation from ideal mixing behavior, which can also be rationalized in terms of electrostatics. Finally, we use a variety of models to investigate the thermal conductivity in the alloys, and we show that the strain field introduced by inhomogeneous chemical bonding is likely to have a significant impact on the transport. The detailed microscopic insight from these studies highlights the benefits of this theoretical modelling approach for studying alloy systems, providing the potential for its extension to a wider range of solid-state materials for thermoelectric applications.