Over the past few decades, seismic studies have revealed complex structural anomalies in the Earth’s deep interior at various scales, such as large low-shear-velocity provinces (LLSVPs) and ultra-low velocity zones (ULVZs) in the lowermost mantle, and small-scale scatterers in the mid-mantle. These structures which are critical for better understanding of the geodynamics and evolution of the deep Earth, need to be further resolved by high-resolution imaging techniques. The spectral-element method (SEM) can be used to accurately simulate seismic wave propagation in heterogeneous Earth models, and its application in full-waveform inversion (FWI) provides a promising high-resolution and high-fidelity imaging technique. But it can be computationally prohibitive when used to model small scale structures in the deep Earth based upon high-frequency seismic waves. The heavy computational cost can be circumvented by using hybrid methods, which restrict the main computation by SEM solver to only a small target region (e.g. above the CMB) encompassing possible 2-D/3-D anomalies, and apply efficient analytical or numerical methods to calculate the wavefield for 1-D background models. These forward modelling tools based on hybrid methods can be then used in the so-called ‘box tomography’ approach to resolve fine-structures in the deep Earth. In this study, we outline the theory of a hybrid method used to model small scale structures in the deep Earth and present its implementation based on SEM solvers in a three-step workflow. First, the wavefield generated by the source is computed for the 1-D background model with traction and velocity saved for the virtual boundary of the target region, which are then used as boundary inputs to simulate the wavefield in the target region based on absorbing boundary condition in SEM. In the final step, the total wavefield at receivers is reconstructed based upon the total wavefield on the virtual boundary computed in the previous step. As a proof-of-concept study, we demonstrate the workflow of the hybrid method based on a 2-D SEM solver. Examples of the hybrid method applied to a coupled fluid–solid model show that our workflow can accurately recover the scattered waves back to the surface. Furthermore, we benchmark the hybrid method on a realistic heterogeneous Earth model built from AK135-F and show how teleseismic scattered waves can be used to model deep Earth structures. By documenting the theory and SEM implementation of the hybrid method, our study lays the foundation for future two-way coupling of 3-D SEM solver with other efficient analytic or numerical 1-D solvers.