Owing largely to multiscale heterogeneity in the underlying fibrous structure, the physics of fluid flow in and through fibrous media is incredibly complex. Using fully resolved finite element (FE) simulations of Newtonian, incompressible fluid flow perpendicular to the fibers, the macroscopic permeability is calculated in the creeping flow regime for arrays of random, ideal, perfectly parallel fibers. On the micro-scale, several order parameters, based on Voronoi and Delaunay tessellations, are introduced to characterize the structure of the randomly distributed, parallel, non-overlapping fiber arrays. In particular, by analyzing the mean and the distribution of the topological and metrical properties of Voronoi polygons, we observe a smooth transition from disorder to (partial) order with decreasing porosity, i.e., increasing packing fraction. On the macro-scale, the effect of fiber arrangement and local crystalline regions on the macroscopic permeability is discussed. For both permeability and local bond orientation order parameter, the deviation from a fully random configuration can be well represented by an exponential term as function of the mean gap width, which links the macro- and the micro-scales. Finally, we verify the validity of the, originally, macroscopic Darcy's law at various smaller length scales, using local Voronoi/Delaunay cells as well as uniform square cells, for a wide range of porosities. At various cell sizes, the average value and probability distributions of macroscopic quantities, such as superficial fluid velocity, pressure gradient or permeability, are obtained. These values are compared with the macroscopic permeability in Darcy's law, as the basis for a hierarchical upscaling methodology.