In this study, 3D finite element analyses (FEA) are conducted to quantify the orthotropic elastic properties and investigate the load transfer mechanism of bone at the sub-lamellar level. Three finite element (FE) unit cells with periodic boundary conditions are presented to model a two-scale microstructure of bone including a mineralized collagen fibril (MCF), the extrafibrillar matrix (EFM) and the resulting fibril array (FAY) under arbitrary loading. The axial and transverse elastic properties of the FAY computed by FEA are calibrated with unique experimental results on ovine micro-samples showing a coherent fibril orientation. They are then systematically compared with those calculated using analytical methods including the basic Voigt, Reuss and shear-lag models, the Mori–Tanaka scheme and the upper and lower bounds by Hashin and Shtrikman. The predicted axial strain ratios between the two-scales are discussed with respect to a recent small-angle X-ray scattering and wide-angle X-ray diffraction study. Beyond apparent elastic properties, the FE models provide stress distributions at both hierarchical levels, confirm the shear lag mechanisms within the MCF and between MCF and EFM and identify potential damage sites under arbitrary loading conditions. A comprehensive sensitivity analysis shows that mineral volume fraction in the fibril array is the dominant parameter on the axial and transverse elastic moduli, while the MCF volume fraction in FAY is the most sensitive variable for the ratio of axial versus transverse elastic modulus followed by the elastic moduli of hydroxyapatite and collagen. The FE model of the FAY developed and calibrated in the current study represents an anatomically realistic, experimentally validated and computationally efficient basis for investigating the apparent yield, post-yield and failure behaviors of lamellar bone in future research.