This study investigates the effect of fracture lower surface roughness on the nonlinear flow behaviors of fluids through fractures when the aperture fields are fixed. The flow is modeled with hydraulic pressure drop = 10 − 4 ~ 10 5 Pa / m by solving the Navier-Stokes equations based on rough fracture models with lower surface roughness varying from JRC = 1 to JRC = 19 . Here, JRC represents joint roughness coefficient. The results show that the proposed numerical method is valid by comparisons between numerically calculated results with theoretical values of three parallel-plate models. With the increment of hydraulic pressure drop from 10-4 to 105 Pa/m spanning ten orders of magnitude, the flow rate increases with an increasing rate. The nonlinear relationships between flow rate and hydraulic pressure drop follow Forchheimer’s law. With increasing the JRC of lower surfaces from 1 to 19, the linear Forchheimer coefficient decreases, whereas the nonlinear Forchheimer coefficient increases, both following exponential functions. However, the nonlinear Forchheimer coefficient is approximately three orders of magnitude larger than the linear Forchheimer coefficient. With the increase in Reynolds number, the normalized transmissivity changes from constant values to decreasing values, indicating that fluid flow transits from linear flow regimes to nonlinear flow regimes. The critical Reynolds number that quantifies the onset of nonlinear fluid flow ranges from 21.79 to 185.19.