The advantage of a smooth representation of bathymetry in terrain-following σ-coordinate ocean models is compromised by the need to avoid numerical errors on steep slopes associated with pressure gradient discretization or spurious diapycnal diffusion. Geopotential z-coordinate models avoid these errors, but greatly underrepresent the interaction of flow with a topographic slope, especially when the bathymetry is underresolved. Hybrid coordinate models are also deficient because it is difficult to find a satisfactory compromise between z and σ coordinates. With volume penalization, we do not seek a compromise, but rather a correction to the usual coordinate systems that realistically recovers continuous and steep bathymetry. The Brinkman volume penalization method studied here is a modified version of the one introduced in Debreu et al. (2020) that simplifies the numerical implementation of the penalization, increases robustness and improves its computational performance for realistic long-term simulations, while preserving accuracy. We apply this penalization method to the Gulf Stream separation problem that has puzzled modelers for decades. The method improves the representation of the flow-topography interaction and achieves realistic separation of the Gulf Stream at resolutions as coarse as 1/8◦. In addition, it provides a tool to separate the effect of eddy activity and topographic slope when changing grid resolution. This has never before been possible because at coarse resolution none of the usual coordinate systems can properly represent a steep continental slope. Our results show that realistic bathymetry is more important than eddy activity in ensuring realistic Gulf Stream separation, even though many recent studies tend to focus on the eddy activity. A steep slope can exert a stabilizing influence that promotes a strong mean slope current with strong inertia that helps it separate from the coast at the topographic curvature of Cape Hatteras. We anticipate that a successful topographic slope correction will be very valuable to climate models, as their current resolution is far from sufficient to represent western boundary currents (WBCs) using traditional coordinate systems. Our results suggest that a climate model with a 1/4◦resolution using volume penalization — and perhaps also some parameterization of the eddy-mean flow interaction to energize the WBCs — would represent ocean circulation much more realistically.