Abstract: Aiming to solve the nonlinear flow in fractured porous media, the coupling characteristics between Darcy flow in pores and Forchheimer flow in fractures are described by means of the pressure transfer function. The finite volume numerical form of seepage equations is derived, and the corresponding numerical code is written. The flow solution by the proposed method for single fracture and intersecting fracture is verified against Frih and Arraras’ solution. Based on this method, the fluid flow behavior of a fractured rock deep-buried tunnel is simulated, which shows it has strong applicability to flow in complex fracture system. The nonlinear flow of tunnel is also analyzed. The results show that the hydraulic gradient of surrounding rock is characterized by "large at bottom and small at top", with the maximum difference of 2.5 times. Therefore, the flow rate at the bottom of the tunnel is greater than that at the top. The distribution homogeneity and density of fracture are the important factors that affect the hydraulic behavior of fractured rock tunnels. At certain water pressure, the more fractures concentrated in the direction of water pressure and the greater the density is, the greater the surrounding rock conductivity is and the greater the flow rate of tunnel is. In this condition, water-inflow accident of tunnels will be prone to occur. The research results may provide reference for the waterproof design and engineering practice of fractured rock tunnels.