Most methods to interpret magnetotelluric data assume that the source fields are plane electromagnetic waves normally incident upon a flat earth. Previous investigations have demonstrated that rough interfaces, such as topography and irregularly layered substructures, can significantly modify magnetotelluric data. In the present study, a modified version of the original Rayleigh scattering theory (Rayleigh-FFT) is extended to an unlimited number of interfaces. With this extension, the Rayleigh-FFT technique allows detailed earth models (including topography) to be used as input for simulation of magnetotelluric fields. This method implicitly includes the contributions of the scattered fields as well as the effects of galvanic and inductive phenomena. Theoretical formulations for three-dimensional structures are presented, but the numerical computations are for two-dimensional magnetotelluric modeling. The uniqueness of using this approach for modeling is that the input is represented by continuous interface points and their local slopes. Therefore, earth models with smoothly varying slopes are effectively calculated. This is in contrast to finite element or difference methods where the earth models are constructed by using rectangular or triangular gridded meshes. There is an inherent error in the Rayleigh representation of the scattered fields which limits the maximum slope of the interface. However, for multiple interfaces, the inherent error is formally eliminated within each rough layer. This error reduction is evidenced by valid solutions at significantly higher slope values relative to previous single interface investigations. Comparisons of the Rayleigh-FFT technique with other independent modeling methods are in agreement where they are considered valid. Independent self-consistency checks are developed to determine the validity of the Rayleigh-FFT solutions. The root-mean- square difference between successive Rayleigh-FFT calculations of apparent resistivity along the surface appears to be the most sensitive diagnostic error indicator. A 1 Ohm-m cosine bell shaped conductor of variable upper interface slope in 100 Ohm-m host rock was investigated for 0.001 Hz. For this model, valid solutions were obtained for slopes up to 1.9 (~62°) and 1.6 (~58°) for transverse electric and transverse magnetic polarizations, respectively. There are still many areas for further investigations concerning the capabilities and limitations using the Rayleigh-FFT technique. However, with the extension to an unlimited number of interfaces, this technique appears to be very promising for magnetotelluric modeling.