This article presents the practical theory that was used to implement the Zeeman effect using Stokes formalism in the Atmospheric Radiative Transfer Simulator ARTS. ARTS now treats the Zeeman effect in a general manner for several gas species for all polarizations and takes into account variations in both magnetic and atmospheric fields along a full 3D geometry. We present how Zeeman splitting affects polarization in radiative transfer simulations and find that the effect may be large in Earth settings for polarized receivers in limb observing geometry. We find that not taking a spatially varying magnetic field into account can result in absolute errors in the measurement vector of at least 10 K in Earth magnetic field settings. The article also presents qualitative tests for O2 lines against previous models (61.15 GHz line) and satellite data from Odin-SMR (487.25 GHz line), and the overall consistency between previous models, satellite data, and the new ARTS Zeeman module seems encouraging.