We address the large-scale regularity theory for the stationary Navier-Stokes equations in highly oscillating bumpy John domains. These domains are very rough, possibly with fractals or cusps, at the microscopic scale, but are amenable to the mathematical analysis of the Navier-Stokes equations. We prove a large-scale Calderón-Zygmund estimate, a large-scale Lipschitz estimate, and large-scale higher-order regularity estimates, namely, C 1,γ and C 2,γ estimates. These nice regularity results are inherited only at mesoscopic scales, and clearly fail in general at the microscopic scales. We emphasize that the large-scale C 1,γ regularity is obtained by using first-order boundary layers constructed via a new argument. The large-scale C 2,γ regularity relies on the construction of second-order boundary layers, which allows for certain boundary data with linear growth at spatial infinity. To the best of our knowledge, our work is the first to carry out such an analysis. In the wake of many works in quantitative homogenization, our results strongly advocate in favor of considering the boundary regularity of the solutions to fluid equations as a multiscale problem, with improved regularity at or above a certain scale. 1. Introduction 171 2. Estimates for the nonlinearity 181 3. Large-scale Lipschitz estimate 193 4. Boundary layers in bumpy John domains 203 5. Higher-order regularity 215 Appendix A. Bogovskii's lemma and some applications 225 Appendix B. Large-scale estimates for the Green's function 227 Appendix C. Proof of the iteration lemma 237 Acknowledgements 239 References 239 MSC2020: 35Q30,