The microstructure of a highly laminated Lower Bowland Shale sample is characterized at the micron-to millimeter scale, to investigate how such characterization can be utilized for microstructure-based modelling of the shale's geomechanical behavior. A mosaic of scanning electron microscope (SEM) back-scattered electron (BSE) images was studied. Mineral and organic content and their anisotropy vary between laminae, with a high variability in fracturing and multi-micrometer aggregates of feldspars, carbonates, quartz and organics. The different microstructural interface types and heterogeneities were located and quantified, demonstrating the microstructural complexity of the Bowland Shale, and defining possible pathways for fracture propagation. A combination of counting-box, dispersion, covariance and 2D mapping approaches were used to determine that the total surface of each lamina is 3 to 11 times larger than the scale of heterogeneities relative to mineral proportion and size. The dispersion approach seems to be the preferential technique for determining the representative elementary area (REA) of phase area fraction for these highly heterogeneous large samples, supported by 2D quantitative mapping of the same parameter. Representative microstructural models were developed using Voronoï tessellation using these characteristic scales. These models encapsulate the microstructural features required to simulate fluid flow through these porous Bowland Shales at the mesoscale.