This paper formulates two different boundary-element+Born series schemes for wave propagation simulation in multilayered media by incorporating a Born series and boundary integral equations.The first scheme directly decomposes the resulting boundary integral equation matrix into the self-interaction operators associated with each boundary itself and the extrapolation operators expressing cross-interactions between different boundaries in a subregion.For the second scheme,the matrix dimension is firstly reduced to a half by the elimination of the traction field in the equations.The resulting new matrix can also be split into the self-interaction matrices associated each subregion itself and the extrapolation matrices interpreting cross-interactions between different subregions in a whole model.Both the numerical schemes avoid the inversion of the relatively much larger boundary integral equation matrix of a full-waveform BE method and hence save computing time and memory greatly.The two schemes are validated by calculating synthetic seismograms for a homogeneous layered model,compared with the full-waveform BE numerical solution.Numerical experiments indicate that both the BEM+Born series modeling schemes are valid and effective.The tests also confirm that the second modeling scheme has a faster convergence in comparison with the first one. This paper formulates two different boundary-element + Born series schemes for wave propagation simulation in multilayered media by incorporating a Born series and boundary integral equations. First solution directly decomposes the resulting boundary integral equation matrix into the self-interaction operators associated with each boundary itself and the extrapolation operators expressing cross-interactions between different boundaries in a subregion. For the second scheme, the matrix dimension is first reduced to a half by the elimination of the traction field in the equations. the new matrix can also be split into the self-interaction matrices associated each subregion itself and the extrapolation matrices interpreting cross-interactions between different subregions in a whole model. But the numerical schemes avoid the inversion of the relatively much larger boundary integral equation matrix of a full-waveform BE method and hence save computing time and memory greatly.The two schemes are validated by calculating synthetic seismograms for a homogeneous layered model, compared with the full-waveform BE numerical solution. Numerical experiments indicate that both both BEM + Born series modeling schemes are valid and effective. The tests also confirm that the second modeling scheme has a faster convergence in comparison with the first one.