A hybrid operator splitting method is developed for computations of twodimensional transverse magnetic Maxwell equations in media with multiple random interfaces. By projecting the solutions into the random space using the polynomial chaos (PC) projection method, the deterministic and random parts of the solutions are solved separately. The novelty of this article lies in using level set functions to characterize the random interfaces and, under reasonable assumptions on the random interfaces, the dimensionality issue from the PC expansions is resolved. For the computational accuracy and efficiency, we will adopt the time-dependent polynomial chaos methods which will be mentioned briefly in this talk.