Inversion of electromagnetic induction data is often realized using smoothing operators for regularizing the process. These operators are applied using derivative matrices, penalizing differences between neighboring mesh cells. The resulting smoothing effect is similar to that of low-pass Gaussian filters and may cause difficulties due to the smoothed boundaries between structures. To prevent such ambiguities, non-linear rank filters are applied widely in image processing applications to filter low frequency variations while keeping boundaries. Correspondingly, we have defined a non-linear rank order smoothing operator for regularization. The effects of the defined parameters are demonstrated on synthetic models and the operator is observed to be able to provide better defined boundaries. The implementation of the operator on field data is realized using single-frequency VLF-R data collected on profiles traversing the surface rupture of the North Anatolian Fault in Kocaeli, Turkey. In the study area, ERT data are also collected on a single profile, and its smooth inversions are realized individually and jointly with VLF-R data for comparisons. Inversions show that the recovered subsurface structure using the rank order smoothing is closer to that indicated by the ERT data and previous studies. According to the recovered models, the clayey fillings in the fault rupture are acting as a barrier affecting the horizontal fluid flow in the near-surface and the resulting water accumulations around the fault rupture may cause susceptibility for earthquake induced damages.