Complex geologic structures generally consist of irregular subregions with piecewise constant properties. Two different boundary-element method (BEM) plus Born-series schemes have been formulated for wave-propagation simulation in such piecewise homogeneous media by incorporating a Born series and boundary integral equations. Both schemes decompose the resulting boundary integral equation matrix into two parts: (1) the self-interaction operator, handled with a fully implicit BEM, and (2) the extrapolation operator, approximated by a Born series. The first scheme associates the self-interaction operator with each boundary itself and interprets the extrapolation operator as cross-interactions between different boundaries in a subregion. The second scheme relates the self-interaction operator toeach subregion itself and expresses the extrapolation operator as cross-interactions between different subregions in a whole model. In the second scheme, the matrix dimension of the resulting boundary integral equation is reduced by eliminating the traction field. Both numerical schemes have been validated by dimensionless frequency responses to a semicircular alluvial valley, compared with the full-waveform BEM numerical solution. We then extended these schemes to a complex fault model by calculating synthetic seismograms to evaluate approximation accuracies. Numerical experiments indicate that the series modeling schemes significantly improve computational efficiency, especially for high frequencies and with multisource seismic surveys. The tests also confirmed that the second modeling scheme has a faster convergence but may cause more noise in higher-order iterations than the first scheme.