A surface integral representation of Maxwell's equations allows the efficient electromagnetic (EM) modeling of three-dimensional structures with a two-dimensional discretization, via the boundary element method (BEM). However, existing BEM formulations either lead to a poorly conditioned system matrix for multiscale problems, or are computationally expensive for objects embedded in layered substrates. This article presents a new BEM formulation which leverages the surface equivalence principle and Buffa-Christiansen basis functions defined on a dual mesh, to obtain a well-conditioned system matrix suitable for multiscale EM modeling. Unlike existing methods involving dual meshes, the proposed formulation avoids the double-layer potential operator for the surrounding medium, which may be a stratified substrate requiring the use of an advanced Green's function. This feature greatly alleviates the computational expense associated with the use of Buffa-Christiansen functions. Numerical examples drawn from several applications, including remote sensing, chip-level EM analysis, and metasurface modeling, demonstrate speed-ups ranging from 3x to 7x compared to state-of-the-art formulations.