Using three-dimensional (3-D) distinct element modeling, we explored a variety of simulations to characterize and interpret the stress permutations in brittle tectonics. Stress inversions of fault slip data or earthquake focal mechanisms often revealed such permutations. The main aim of our study is to produce simple, mechanically consistent 3-D models that account for these switches between the principal stress axes σ1/σ2 or σ2/σ3. Even with simple boundary conditions the stress changes induced by variations in rheology are large enough to modify the local tectonic behavior and to produce permutations of principal stress axes. Rather than simple directional changes of stress axes, which exist but often remain limited, the relative variations in principal stress magnitudes are the major cause of permutations σ1/σ2 and σ2/σ3. In nature, permutations being left apart, the orientations of axes often remain tightly clustered. In our experiments we adopted a ratio Φ = (σ2 - σ3)/(σ1 - σ3) of 0.5, which makes permutations difficult (low Φ favors σ2/σ3 permutations, high Φ favors σ1/σ2 permutations), and we explored a variety of tectonic situations involving compression, extension, and strike slip. Our experiments indicate that the major causes of stress permutations are the heterogeneity of the brittle deformation (e.g., intact rock massifs between heavily faulted deformation zones) and the anisotropy of the mechanical properties that results from fracturing and faulting, which concur to modify the mechanical balance inside the analyzed volume and to produce stress permutations. Our models show that contrasts and anisotropy in rock properties favor stress permutations. Of major importance is the existence of relatively resistant zones at tips of the deformed zones, acting as channels where stress concentrates and switches occur. Because in nature such zones move in time and space, it is not surprising that stress permutations are so pervasive.