A new surface integral formulation designed to efficiently simulate effects of 3D topography and displacement currents on radiomagnetotelluric data is presented. Surface charges, magnetic surface currents and electrical surface currents are introduced to ensure the stability over a wide rage of frequencies. Using a triangular parameterization and the Galerkin approach, the surface integral formulation is transformed into a linear system of equations, in which several double surface integrals must be evaluated. If the triangles in the outer integrals are close to, adjoin or overlap with the triangles in the inner integrals, then these double surface integrals will be weakly or strongly singular. An accurate and robust singularity extraction technique is applied to transform these singular integrals, such that they can be evaluated by analytical expressions. The remaining logarithmic singularities in the outer integrals are removed with a smart strategy of interchanging the order of integration. If the triangles fully overlap, strongly singular integrals can be transformed into weakly singular ones using Gauss' divergence theorem. The efficacy of our singularity extraction technique is validated by means of numerical examples.