Hybrid finite element (FE)-boundary integral (BI) analysis of infinite periodic arrays is extended to include planar multilayered Green's functions. In this manner, a portion of the volumetric dielectric region is modeled via the FE method whereas uniform multilayered regions are characterized using a multilayered Green's function. As such, thick uniform substrates can be modeled without loss of efficiency and accuracy. The multilayered Green's function is analytically computed in the spectral domain and the resulting BI matrix-vector products are evaluated via the fast spectral domain algorithm (FSDA) without explicit generation of the BI-matrix. Furthermore, the number of Floquet modes in the BI expansion is kept very few by appropriately placing the BI surfaces within the computational unit cell. As a result, the computational cost of the method is very little and the model can easily be adapted to various requirements. Examples of frequency selective surface (FSS) arrays are analyzed with this method to demonstrate the accuracy and capability of the approach.