Seismic tomography has long been recognized as one of the most effective and widely applied tools for probing the Earth’s interior. Since the pioneering studies of K Aki and collaborators in the 1970s, advances in seismic observational techniques, computational capabilities, and forward and inversion methods have led to remarkable images of the Earth’s internal structure, greatly improving our understanding of the Earth’s compositional status and dynamic process. Ray-theoretical traveltime tomography successfully revealed the large-scale features of the lithosphere and mantle. However, the high-frequency approximation inherent in traveltime tomography prevented the full usage of information embedded in seismic records and thus limited the achievable resolution with a given station coverage. This, combined with the advances in computational capabilities, naturally gave rise to the development of finite-frequency approach and ultimately the full-waveform paradigm so that all structural information embedded in seismic records can be fully employed to invert any model parameters such as wave velocity, density, anisotropy, and anelasticity. A recent development employs full waveforms of short-period teleseismic arrivals recorded by a local seismic array to image the subarray structures. Since the approach only resolves the structures within the volume beneath the array, it is also referred to as “box tomography”. It achieves both high computational efficiency and spatial resolution, and has shown to be effective in resolving lithospheric structures when applied to dense arrays, particularly in tectonically complex but seismically inactive regions. Moreover, it offers robust constraints on multiple structural parameters, including density, P- and S-wave velocities, and even anisotropy.

In box tomography, the teleseismic wavefields are calculated by a hybrid approach, in which the wavefields outside the box of interest are calculated by efficient algorithms using simplified （often 1-D） models, and this external wavefields are injected into the subarray box for wavefield simulation by accurate wave equation solvers in complex 3-D models that incorporates realistic topography and heterogeneities. The sensitivity kernels of the full-waveform objective function to model perturbations are computed through the adjoint method that involves a forward simulation followed by an adjoint simulation, in which the data residuals are back-propagated from the stations. The resulting forward and adjoint wavefields are convolved to construct the gradients of the objective function with respect to the model parameters such as density and P- and S-wave velocities. The inversion is conducted iteratively using an efficient optimization scheme.

Practical applications of box tomography generally make use of teleseismic records from densely spaced seismic arrays. The teleseismic events with high waveform quality are selected and preprocessed through component rotation, normalization, windowing, and the hierarchical multi-frequency filtering to reduce the non-linearity of the full-waveform objective function. An apparent source time function is estimated for each teleseismic event by deconvolving the synthetic seismograms from the corresponding records, thereby mitigating the influence of the inaccuracies in both the structural model outside the target box and the earthquake source model. The spatial resolution of the resulting model inside the target box can be evaluated through resolution tests such as recovery test, checkerboard test and spike test, which provide illustrative indications of the model reliability. These procedures have been applied and proven effective in several box tomographic studies on orogenic belts, subduction zones, and intra-continental crust and upper-mantle structures, demonstrating the capability of box tomography to resolve small-scale lateral heterogeneities in the lithosphere in diverse tectonic environments.

The known limitations and challenges of the box tomography method include: ① the occurrence of singularities in the sensitivity kernels near seismic stations, which arises from the analytic singularities of the adjoint wavefields at the receivers; ② the inaccuracies in the global model outside the box of interest used for computing the teleseismic wavefields to propagate into the box, leading to errors in the hybrid simulation results; ③ potential influence of errors in the teleseismic earthquake source model, which is assumed known and kept fixed during the inversion.

Future directions of box tomography development revolve around the three limitations and challenges. First, it is necessary to design effective ways to suppress the singularities near stations. Possible solutions include the introduction of station-based weighting and data-adaptive model discretization. Second, using more accurate （2-D or even 3-D） Earth models for calculating the teleseismic wavefields in hybrid simulations. Third, the source-related errors may be mitigated by adopting objective functions less dependent on earthquake sources, for instance, the receiver function. In addition, joint inversion of teleseismic body and surface waves should be further explored. Such approaches exploit the complementary sensitivity of different wave types, yielding improved depth resolution in imaging crustal and upper-mantle structures.

In this paper, we present a comprehensive review on box tomography. Our aim is to spark interests and facilitate applications of box tomography, which serves as an effective tool for the investigation of the earth’s lithospheric structure.