The distortion of atomic structure around In and Cu dopants in sphalerite ZnS was explored by extended X-ray absorption fine structure (EXAFS) spectroscopy enhanced by reverse Monte Carlo (RMC) simulation (RMC-EXAFS method). These data were complemented with quantum chemical Density Functional Theory (DFT) calculations and theoretical modeling of X-ray absorption near edge spectroscopy (XANES) spectra. The RMC-EXAFS method showed that in the absence of Cu, the In-bearing solid solution is formed via the charge compensation scheme 3Zn2+↔2In3+ + □, where □ is a Zn vacancy. The coordination spheres of In remain undistorted. Formation of the solid solution in the case of (In, Cu)-bearing sphalerites follows the charge compensation scheme 2Zn2+↔Cu+ + In3+. In the solid solution, splitting of the interatomic distances in the 2nd and 3rd coordination spheres of In and Cu is observed. The dopants’ local atomic structure is slightly distorted around In but highly distorted around Cu. The DFT calculations showed that the geometries with close arrangement (clustering) of the impurities—In and Cu atoms, or the In atom and a vacancy—are energetically more favorable than the random distribution of the defects. However, as no heavy In atoms were detected in the 2nd shell of Cu by means of EXAFS, and the 2nd shell of In was only slightly distorted, we conclude that the defects are distributed randomly (or at least, not close to each other). The disagreement of the RMC-EXAFS fittings with the results of the DFT calculations, according to which the closest arrangement of dopants is the most stable configuration, can be explained by the presence of other defects of the sphalerite crystal lattice, which were not considered in the DFT calculations.