Viscoplastic non-Newtonian fluid flow in rock fractures is critical for subsurface applications but remains understudied. In this work, the flow of Bingham plastic fluids in a single rough fracture is simulated and analyzed in detail. For this purpose, the yield stress rheology of kaolinite suspensions at concentrations of 20 and 28.5 wt. % is modeled using the Bingham-Papanastasiou framework. Simulations are carried out over a broad range of incoming velocities (i.e., u( 0 )= 0.001-2 m s(-1)), covering flow regimes from negligible to strong inertia effects. The results reveal that the yield stress effect creates rigid zones, either moving with the flow or attached to the fracture surface, while increasing the effective viscosity and, therefore, stabilizing the fracture flow. At high flow velocities, this effect diminishes, reducing the area of rigid zones and promoting nonlinear features such as the manifestation of multi-scale eddies. Flow tortuosity increases with the incoming velocity but decreases with yield stress, reflecting the stabilizing influence of viscoplasticity. The total pressure drop demonstrates non-Darcy behavior for both Newtonian and yield stress fluids, driven by the combined effects of fracture roughness and yield stress. Forchheimer's equation effectively predicts the pressure drop across the fracture, capturing nonlinear contributions.