ABSTRACTThe description of chemical phenomena in solution is as challenging as it is important for the accurate calculation of molecular properties. Here, we present the implementation of the polarisable continuum model (PCM) in the four-component Dirac–Kohn–Sham density functional theory framework, offering a cost-effective way to concurrently model solvent and relativistic effects. The implementation is based on the matrix representation of the Dirac–Coulomb Hamiltonian in the basis of restricted kinetically balanced Gaussian-type functions, exploiting a non-collinear Kramer's unrestricted formalism implemented in the program ReSpect, and the integral equation formalism of the PCM available through the stand-alone library PCMSolver. Calculations of electron paramagnetic resonance parameters (g-tensors and hyperfine coupling A-tensors), as well as of the temperature-dependent contribution to paramagnetic nuclear magnetic resonance (pNMR) shifts, are presented to validate the model and to demonstrate the importance of taking both relativistic and solvent effects into account for magnetic properties. As shown for selected Ru and Os complexes, the solvent shifts may amount to as much as 25% of the gas-phase values for g-tensor components and even more for pNMR shifts in some extreme cases.