The anisotropy of surface roughness in shale fractures affects gas flow characteristics, which are crucial to evaluating shale gas production. However, the coupled relationship among surface roughness distribution, effective stress, and gas flow in a shale fracture has not been quantified sufficiently. This study numerically investigated the nonlinear gas flow in shale anisotropic fracture under varying effective stresses. The results show that the anisotropy of the roughness characteristics for shale fracture is characterized by autocorrelation function. The stress field distribution and flow characteristics on the fracture show distinct layering patterns: hourglass layering for stress, parallel layering for pressure, and elliptical layering for velocity. As shale gas transports from profile lines of x = 1 mm to x = 85 mm, the average pressure loss rate increment increases from 0.055% MPa −1 to 3.875% MPa −1 , demonstrating an enhanced effect of effective stress on the flow with increasing x value. Furthermore, in microscale shale fracture, flow velocity variation dominates in determining flow rate compared to fracture aperture variation during loading. Finally, as effective stress increases, fracture closure variation among different profile lines becomes increasingly significant, while relative differences in friction coefficient and Reynolds number among lines show a negative correlation with effective stress.