We characterize seepage flow in porous media via a series of three-dimensional (3D) direct numerical simulations (DNSs) of Navier-Stokes flows in representative volume elements (RVEs) at the pore level. The immersed boundary method with a fixed staggered grid is employed for calculations based on the finite difference method in the pore domain formed with spatially arranged rigid particles in RVEs. The numerical results of the DNSs with different particle sizes and different seepage flow velocities are volume-averaged over the RVEs to evaluate the permeability coefficients of the seepage flows and are ordered according to the Reynolds numbers for porous media. After the numerical scheme is validated with a simple RVE model, the transition from Darcy to non-Darcy flow, for which experimental results have reported in the literature, is demonstrated in terms of the calculated permeability coefficients and examined based on the pore-level flow characteristics to clarify the mechanism of the permeability change, i.e. the macroscopic behaviour of RVEs related to the Reynolds number.