A methodology is presented for the Krylov subspace-based model order reduction of finite element models of electromagnetic structures with material properties and impedance boundary conditions exhibiting arbitrary frequency dependence. The proposed methodology is a generalization of an equation-preserving Krylov model order reduction scheme for methodology for second-order, linear dynamical systems. The emphasis of this paper is on the application of this method to the broadband model order reduction of planar circuits including lossy strips of arbitrary thickness and lossy reference planes. In particular, it is shown that the proposed model order reduction methodology provides for the accurately modelling of the impact of the frequency dependence of the internal impedance per unit length of the thick lossy strip on the electromagnetic response of the stripline structure over a very broad, multi-GHz frequency band, extending all the way down to frequencies in the DC neighbourhood. In addition, the application of the proposed methodology to the broadband modelling of electromagnetic coupling between strips on either side of a lossy ground plane is demonstrated.