Chao Lyu

Liang Zhao
email: zhaoliang@mail.iggcas.ac.cn

Yann Capdeville

Novel hybrid numerical simulation of the wave equation by combining physical and numerical representation theorems and a review of hybrid methodologies

Hybrid wave numerical simulation methods of circumventing the heavy compu-10 tational cost in the global waveform tomography are reviewed.

11

• The proposed hybrid method has the flexible local meshing and is highly accu-12 rate and memory efficient.

13

• Only three physical quantities located exactly on the hybrid interface are required 14 to construct the hybrid inputs.

Introduction

Over the past three decades, with the development of effective computer clusters, efficient numerical methods, such as the finite difference method (FDM) and spectral element method (SEM), improved full waveform inversion (FWI) techniques, and exten-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth sive global deployment of stations, the imaging of Earth has achieved unprecedented accuracy [START_REF] French | Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography[END_REF][START_REF] Bozdag | Global adjoint tomography: first-generation model[END_REF][START_REF] Fichtner | The collaborative seismic earth model: Generation 1[END_REF][START_REF] Lei | Global adjoint tomography-model GLAD-M25[END_REF]. Faster and more accurate calculations of the elastic wavefield within heterogeneous media are key to developing and improving accurate imaging techniques relying on full-waveform analysis [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Pratt | Gaussnewton and full newton methods in frequencyspace seismic waveform inversion[END_REF][START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF][START_REF] Tromp | Seismic wavefield imaging of earths interior across scales[END_REF][START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF]. Earths structure is multi-scale, yet capturing such a broad range of complex heterogeneities with seismic wave propagation across the observable frequency band (e.g., ≥ 1 s) requires thousands of global numerical simulations of seismic/acoustic wave equations, which are still computationally prohibitive. Therefore, it is important to lighten this computational burden to further increase the resolution and accuracy of imaging. Many research efforts have been devoted to speed up global wave simulations. For example, the real coupling "CSEM" method couples the SEM simulation with normal mode calculations [START_REF] Capdeville | Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models[END_REF]. Another two-step approach is based on the physical domain decomposition, where wave propagation is first performed in the global reference model and the hybrid inputs are calculated and saved only once. Subsequently, they are imposed into the local target model as the equivalent virtual sources to perform the hybrid simulation [START_REF] Wen | A two-dimensional P-SV hybrid method and its application to modeling localized structures near the core-mantle boundary[END_REF][START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF][START_REF] Chevrot | Shear wave splitting in threedimensional anisotropic media[END_REF][START_REF] Chen | A wave equation migration method for receiver function imaging: 1. theory[END_REF][START_REF] Zhao | A two-dimensional hybrid method for modeling seismic wave propagation in anisotropic media[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF][START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF]Tong, Chen, et al., 2014;[START_REF] Zhao | Direct numerical modeling of time-reversal acoustic subwavelength focusing[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF][START_REF] Leng | 3-D scattering of elastic waves by small-scale heterogeneities in the Earths mantle[END_REF][START_REF] Pienkowska | High-frequency global wavefields for local 3-D structures by wavefield injection and extrapolation[END_REF]. Hereafter, we focus on such a two-step approach and refer to it as the "hybrid method" (the hybrid simulation between the global reference and local target models). It can be applied for imaging localized multi-scale heterogeneities in the shallow and deep Earth [START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the semdsm hybrid method[END_REF][START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic p waves[END_REF]Masson & Romanowicz, 2017a[START_REF] Masson | Fast computation of synthetic seismograms within a medium containing remote localized perturbations: a numerical solution to the scattering problem[END_REF][START_REF] Zhang | Linear array ambient noise adjoint tomography reveals intense crust-mantle interactions in north china craton[END_REF][START_REF] Beller | Lithospheric architecture of the south-western Alps revealed by multiparameter teleseismic full-waveform inversion[END_REF][START_REF] Clouzet | Box Tomography: first application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[END_REF]K. Wang, Yang, et al., 2021;K. Wang, Wang, et al., 2021), and it is referred to as the "box tomography" (Masson & Romanowicz, 2017a).

Because almost all the forward and backward numerical simulations of wave propagation are confined to the smallest computational volume, local high-frequency numerical simulations can be implemented without difficulty. Moreover, the specially selected spatial illumination from remote seismic events and/or stations facilitates higher-resolution imaging. Its principal difficulty results from unknown mediums in the external domain that mask the target structure. [START_REF] Masson | Fast computation of synthetic seismograms within a medium containing remote localized perturbations: a numerical solution to the scattering problem[END_REF] demonstrated that an unbiased image can be inverted using box tomography if the initial background model is a homogenized true model [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF] in a two-dimensional (2D) global acoustic model framework. These studies provided an important theoretical foundation for the practical applications of box tomography, initially demonstrating its promising

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth application prospects for probing key three-dimensional (3D) structures anywhere within the Earth.

As the forward part of box tomography, an efficient implementation of hybrid numerical simulation is crucial. The first step of the global wave simulation in the hybrid wave simulation framework has gradually evolved from the approximated generalized ray theory (GRT) calculations [START_REF] Wen | A two-dimensional P-SV hybrid method and its application to modeling localized structures near the core-mantle boundary[END_REF][START_REF] Chen | A wave equation migration method for receiver function imaging: 1. theory[END_REF][START_REF] Zhao | A two-dimensional hybrid method for modeling seismic wave propagation in anisotropic media[END_REF][START_REF] Zhao | Upper mantle seismic anisotropy beneath a convergent boundary: SKS waveform modeling in central tibet[END_REF][START_REF] Zhao | First seismic evidence for continental subduction beneath the western alps[END_REF] to the direct solution method (DSM) [START_REF] Kawai | Complete synthetic seismograms up to 2 hz for transversely isotropic spherically symmetric media[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF][START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic p waves[END_REF], planar wave approximation frequencywavenumber (FK) [START_REF] Zhu | A note on the dynamic and static displacements from a point source in multilayered media[END_REF]Tong, Chen, et al., 2014;Tong, Komatitsch, et al., 2014;[START_REF] Monteiller | On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF], AXISEM [START_REF] Nissen-Meyer | Axisem: broadband 3-D seismic wavefields in axisymmetric media[END_REF][START_REF] Beller | Lithospheric architecture of the south-western Alps revealed by multiparameter teleseismic full-waveform inversion[END_REF][START_REF] Pienkowska | High-frequency global wavefields for local 3-D structures by wavefield injection and extrapolation[END_REF][START_REF] Monteiller | On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF], SPECFEM3D globe [START_REF] Komatitsch | The spectral-element method, beowulf computing, and global seismology[END_REF][START_REF] Clouzet | Box Tomography: first application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[END_REF], and AXISEM3D [START_REF] Leng | Axisem3D: broadband seismic wavefields in 3-D global earth models with[END_REF][START_REF] Leng | 3-D scattering of elastic waves by small-scale heterogeneities in the Earths mantle[END_REF], owing to the large computation resources required for 3D high-frequency wave numerical simulations at the global scale.

For the second step of the hybrid wave numerical simulation, a unifying theory was introduced in [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF], the representation theorem being their theoretical basis [START_REF] Wen | A two-dimensional P-SV hybrid method and its application to modeling localized structures near the core-mantle boundary[END_REF][START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF]. Existing hybrid-wave numerical simulation methods consist of two main categories. The first one explicitly approximates the surface integral(s) in the physical representation theorem to build the physical hybrid inputs and is referred to as the "multiple point sources method" [START_REF] Chevrot | Shear wave splitting in threedimensional anisotropic media[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF]Tong, Chen, et al., 2014;Tong, Komatitsch, et al., 2014;[START_REF] Zhao | Direct numerical modeling of time-reversal acoustic subwavelength focusing[END_REF][START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic p waves[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF]. The second category constructs the numerical hybrid inputs using a spatial window function and the discrete wave equation based on the numerical representation theorem and is referred to as the "direct discrete differentiation method" [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF][START_REF] Yoshimura | Domain reduction method for three-dimensional earthquake modeling in localized regions, Part II: verification and applications[END_REF][START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF][START_REF] Masson | Fast computation of synthetic seismograms within a medium containing remote localized perturbations: a numerical solution to the scattering problem[END_REF], 2017a;[START_REF] Clouzet | Box Tomography: first application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[END_REF].

In the research on hybrid wave numerical simulation methods, the following operations are vital: in the first step, reducing the memory requirements of the hybrid inputs in the global simulation; in the second step, accurately performing the hybrid numerical simulations completely independent of the first global simulation (e.g., the irrelevant selections of spatial meshing and temporal steps) and increasing the computation speed of the local hybrid numerical simulations. To the best of our knowledge, no previous hybrid method has combined the two main categories of hybrid simulation. The proposed "new" hybrid method inherits the advantages of both main hybrid methods and simultaneously eliminates their disadvantages. To save hybrid inputs, the memory requirements of the proposed hybrid method and memory-saving multiple point sources method are the same. To use the proposed method in the acoustic scenario, we must know

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth three incident physical quantities: the gradient, potential, and acceleration. Note that the acceleration term can be calculated using the stored potential term during the hybrid simulation; therefore, it does not require storing. It has the advantage of being naturally combined with a perfectly matched layer (PML) for very heterogeneous media without a significant artificial error. The storage ratio between the new hybrid method and the direct discrete differentiation methods based on the numerical representation theorem is approximately 2 N +1 (N : the polynomial degree in the SEM). Therefore, it is more suitable for the hybrid numerical simulations of wave propagation, particularly for an efficient very high polynomial degree SEM [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF].

The remainder of this article is organized as follows. First, we illustrate the detailed expression of the fundamental representation theorem. Subsequently, two main categories of the hybrid simulation, including five different subcategories [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF][START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF][START_REF] Zhao | Direct numerical modeling of time-reversal acoustic subwavelength focusing[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF] are analyzed, and the new hybrid method is proposed based on their combination. Thereafter, a series of 2D and 3D hybrid acoustic wave numerical simulations using SEM in homogeneous and heterogeneous models demonstrate its validity. Finally, the discussion and conclusions are presented.

Methodology

In this section, we begin with the elastic and acoustic wave equations, and we briefly introduce the principles of the spectral element method and the associated representation theorems. Subsequently, some nomenclature and the overall workflow of the hybrid simulation are introduced. Thereafter, we briefly analyze the widely used hybrid methods, including two main categories and five subcategories, and summarize their respective merits and limitations. Finally, the new hybrid method is proposed based on their combination.

Elastic and Acoustic Wave Equation

The propagation of seismic waves is governed by the equations of motion:

ρü = ∇ • σ + f σ = C : ε ε = 1 2 [∇u + (∇u) T ], (1)
where u(x, t) is the displacement field vector, ρ(x) is the density, σ(x) is the stress tensor, ε(x) is the strain tensor, and f (x, t) are the body forces in the elastic domain Ω; u is subject to boundary conditions on ∂Ω (i.e., traction vanishes at the Earth's surface).

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

In the acoustic scenario, considering an acoustic domain Ω, the velocity potential q is the solution of

1 κ q = ∇ • u + f u = 1 ρ ∇q, (2)
where κ(x) is the acoustic bulk modulus, ρ(x) is the (mass) density, u(x, t) is the displacement vector, and f (x, t) is a scalar source term. Generally, the acoustic medium is fully described using only two parameters: the density ρ(x) and the speed of sound

V (x) such that κ(x) = ρ(x)V 2 (x)
. Note that we only tested our proposed hybrid method in the acoustic scenario in the following numerical experiments.

Principle of the Spectral Element Method

The SEM is a finite-element type method. Rather than directly using the boundary conditions and equations of motion, it is based on an associated weak form by dotting the wave equation with an arbitrary test function w(x), integrating by parts over the model domain Ω, and imposing the associated boundary conditions. To numerically solve the weak form of the equation, we require a discrete approximation of the continuous problem. For the spatial discretization of the SEM, the 3D model domain Ω is spatially split into N e non-overlapping spectral elements. For the acoustic wave equation,

q e (x) = i,j,k=0,N +1 Q e ijk w e ijk (x) (e = 1, • • • , N e ), w e ijk (x) is the i + j × (N + 1) + k×(N +1
) 2 th test function in the eth spectral element, q e (x) is q(x) restricted to the eth element, the subscript N represents the polynomial degree of spectral element, and

Q e ijk is q(x e ijk )
where x e ijk is the coordinate of the i+j×(N +1)+k×(N +1) 2 th Gauss-Lobatto-Legendre (GLL) point in the eth element. After the numerical integration based on the GLL points is applied to the weak form of the above acoustic wave equation (2), the following ordinary differential equation (ODE) is obtained:

M Q + KQ = F, (3)
where Q is the discrete potential vector, M is the global diagonal mass matrix, K is the global stiffness matrix, and F is the discrete source vector [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF].

Note that KQ is calculated using the tensor product.

Theoretical Basis of Hybrid Simulation: Representation Theorem

Hybrid methods are often based on the representation theorem [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF][START_REF] Zhao | Direct numerical modeling of time-reversal acoustic subwavelength focusing[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF]. It quantitatively states a manner in which displacement u at a certain point x consists of three contributions because of the force f (x, t) throughout the closed V , as well as contributions owing to

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth the traction T(u, n) and the displacement u itself on the surface S = ∂V .

u n (x, t) = ∞ -∞ dτ V f i (x , τ )G ni (x, t -τ ; x , 0)dV (x ) + ∞ -∞ dτ S T i (u(x , τ ), n)G ni (x, t -τ ; x , 0)dS(x ) - ∞ -∞ dτ S u i (x , τ )C ijkl n j G nk,l (x, t -τ ; x , 0)dS(x ), (4)
where Green's tensor G ni (x, t; x , τ ) denotes the nth direction component of the displacement at (x, t) owing to a unit point force at (x , τ ) in the ith direction,

T i = σ ij n j =
C ijkl u k,l n j is the ith component of the traction on S, n j is the jth component of the normal vector n on S, C ijkl is the fourth-order elasticity tensor, and the notation G nk,l rep-

resents ∂G nk ∂x l
. The third integral in Equation ( 4) indicates the displacement contribution owing to the moment density tensor m kl = u i n j C ijkl . Note that the Einstein summation convention is used. Equation ( 4) indicates that if the exciting force f i (x, t) is known throughout the volume V , the traction T i (x, n, t) and the wavefield u i (x, t) are known on the surface S, then the displacement u n (x, t) within the volume V can be regenerated, which is the theoretical foundation of the hybrid simulation of the elastic wave equation.

Similar to the elastic scenario, the acoustic pressure can be expressed formally in the temporal domain using the Helmoltz-Kirchhoff representation theorem:

q(x, t) = ∞ -∞ dτ V f (x , τ ) G (x, t -τ ; x , 0) dV (x ) + ∞ -∞ dτ S 1 ρ ∇ q (x , τ ) G (x, t -τ ; x , 0) • ndS(x ) - ∞ -∞ dτ S 1 ρ q (x , τ ) ∇ G (x, t -τ ; x , 0) • ndS(x ). (5)
Equation ( 5) can be transferred into the frequency domain as follows (the frequency dependence of functions is omitted hereafter to avoid complexity for equations):

q(x) = V f (x )G(x; x )dV (x ) + S 1 ρ ∇ q(x )G(x; x ) • ndS(x ) - S 1 ρ q(x )∇ G(x; x ) • ndS(x ). ( 6
)
Note that the same notation is used for quantities in the time and frequency domains.

By replacing the surface sources with equivalent body sources [START_REF] Aki | Quantitative seismology[END_REF], Equation (6) can be rewritten as

q(x) = V f (x ) + f T (x ) + f q (x ) G(x; x )dV (x ), (7)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth where the two equivalent body forces associated with the gradient ∇q(x, t) and the potential q(x, t) are given by

f T (x ) = S 1 ρ n i ∂q ∂η i δ(x -η)dS(η) f q (x ) = S 1 ρ ∂ ηi [n i q(η)δ(x -η)]dS(η), (8)
where the terms T = 1 ρ n i ∂q ∂xi and m i = n i q are similar to the traction and moment density tensor terms in the representation theorem for elastic waves. Equation ( 5) indicates that if the exciting force f (x, t) is known throughout the volume V , the gradient ∇q(x, t) and the potential q(x, t) are known on the surface S, then the pressure q(x, t)

within the volume V can be regenerated. Similar to the elastic scenario, it is also the theoretical foundation of the hybrid simulation of the acoustic wave equation.

Nomenclature and Workflow of Hybrid Simulations

For a clearer statement, some nomenclature often used in the hybrid simulations are listed below and depicted in Figure 1.

• Global domain Ω g : the entire physical global domain containing the local domain Ω l and external domain Ω e .

• Local domain Ω l : the local closed box surrounded by the hybrid interface S and located inside the global domain Ω g .

• External domain Ω e : the external part of the global domain Ω g outside the local domain Ω l .

• Global reference model M g0 : the known 1D model (e.g., 1D AK135 [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF]) or 3D long-period velocity model (e.g., SEMUCB-WM1 in [START_REF] French | Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography[END_REF]) from previous studies assigned to the global domain Ω g , which includes the known external model M e0 and local reference model M l0 .

• Global target model M g1 : the global model also assigned to the global domain Ω g , which includes the assumed known external model M e0 and the unknown local target model M l1 .

• Local reference model M l0 : the known model assigned to the local domain Ω l , which is always the starting model for box tomography.

• Local target model M l1 : the unknown model to be inverted for and also assigned to the local domain Ω l .

• External model M e0 : the known external model of the global reference or target model assigned to the external domain Ω e .

• Hybrid interface S: the interface separating the external domain Ω e and local domain Ω l .

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Note that Ω e + Ω l = Ω g for the domains, and M e0 + M l0 = M g0 and M e0 + M l1 = M g1 for the models. As the foundation of box tomography, the forward and backward numerical simulations must be exclusively re-performed in the unknown local target models M l1 because the external model M e is assumed to be unperturbed in the framework of the iterative box tomography.

The key workflow of the hybrid simulation for forward solving the wave equation in box tomography consists of the following three main steps:

• (i) Calculate and record the associated hybrid inputs using global solvers: GRT, FK, DSM, or AxiSEM for the 1D global reference model M g0 (e.g., AK135), and SPECFEM3D globe or AxiSEM3D for the 3D global reference model M g0 (e.g., SEMUCB-WM1).

• (ii) Impose the hybrid inputs recorded in the first step into the local target model M l1 and perform the local hybrid numerical simulation using a local solver, such as the program RegSEM [START_REF] Cupillard | Regsem: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale[END_REF].

• (iii) Record the residual wavefields on the hybrid interface S and use a wavefield extrapolation method [START_REF] Robertsson | An efficient method for calculating finitedifference seismograms after model alterations[END_REF] to obtain the seismic response for the receivers outside the local domain Ω l (Masson & Romanowicz, 2017a[START_REF] Masson | Fast computation of synthetic seismograms within a medium containing remote localized perturbations: a numerical solution to the scattering problem[END_REF].

Multiple Point Sources Method

The first major category, the multiple point sources method, primarily consists of the "VM" (abbreviation obtained from [START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF]) and "RP" (abbreviation for the representation theorem) hybrid methods.

VM Hybrid Method

The first popular type of the multiple point sources method is the VM hybrid method [START_REF] Chevrot | Shear wave splitting in threedimensional anisotropic media[END_REF][START_REF] Godinho | Numerical Simulation of Ground Rotations along 2D Topographical Profiles under the Incidence of Elastic Plane Waves[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF]Tong, Chen, et al., 2014;Tong, Komatitsch, et al., 2014;[START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the semdsm hybrid method[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF]. The concept of the local target model M l1 associated with the VM method is shown in Figure 1c.

The hybrid interface S cannot be selected at will, but it should be the same as the boundary of the local domain ∂Ω l indicted by the closed dashed cyan line. The VM hybrid method is only based on a special representation theorem with the traction-free boundary condition, which will cause the second displacement surface integral term equal to be zero, as discussed on page 29 of the 2nd edition of [START_REF] Aki | Quantitative seismology[END_REF] and Equation ( 12)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth in [START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF]:

q(x, t) = ∞ -∞ dτ S 1 ρ ∇ q (x , τ ) G (x, t -τ ; x , 0) • ndS(x ). (9)
The physical hybrid input is the traction ∇q (gradient) recorded in the global reference model M g0 . Based on the numerical integration of Equation ( 9) as implemented in Equation (A3), the physical hybrid input ∇q can be imposed on the boundary ∂Ω l of the local target model M l1 to implement the hybrid numerical simulation (refer to Appendix A for details).

The VM hybrid method physically operates owing to the effect of the doubling amplitude of the free reflection from the hybrid interface. Note that here, only the first gradient surface integral in Equation ( 6) is imposed on the hybrid interface S, and the second potential integral term is equal to zero. We can consider that a force source is imposed at a point within the local domain, which is infinitely close to the hybrid interface. According to the source mechanism of the point source, the amplitudes of the two wavefronts propagating in opposite directions are the same. One wave (with the halfamplitude of the hybrid waveform) propagating towards the hybrid interface will have an equivalent free reflection after encountering the interface; it then superimposes the other wave (with the other half-amplitude of the hybrid waveform) propagating away from the hybrid interface S and into the local domain. Thus, the free reflection from the hybrid interface has the effect of doubling the amplitude, with full recovery inside the local domain Ω l . This hybrid method does not have an external region and loses the ability to use the PML.

For the hybrid simulation using the VM method, when the local target model M l1 is selected to be same as the local reference model M l0 , the hybrid wavefield (the wavefield of the hybrid numerical simulation in the local domain Ω l ) can be recovered by imposing the physical hybrid inputs, for example, by using the SEM as illustrated in Equation (3):

M 0 Q0 + K 0 Q 0 = F T 0 , (10)
where Q 0 is the reference discrete potential vector, M 0 is the global mass matrix, and K 0 is the global stiffness matrix assigned to the local reference model M l0 . F T 0 is the imposed hybrid source term, which is constructed using Equation (A3) only with the gradient ∇q 0 of some integration points (e.g., GLL points used in this study) exactly on the hybrid interface S. Please note that this is only applicable to the local reference model.

If some perturbations are present in the local domain, then q and ∇q are not actually known on the hybrid interface S because they are different from the global simulation.

Appendix B discusses the related hybrid simulation with a limited absorbing boundary

-10-

Accepted Article

This article is protected by copyright. All rights reserved. 4b2 depict the hybrid wavefield and waveforms calculated using the VM hybrid method in the local reference model M l0 . We observe that the VM hybrid method is sufficiently accurate. The related error was ≈ 0.01% under the spatial meshing setting when the local reference model M l0 was adopted (Section 4.1.1). Note that Figure 4 does not use any absorbing boundary condition outside the hybrid interface.

RP Hybrid Method

The second type of the multiple point sources hybrid method is the RP hybrid method, which is based on the normal representation theorem [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF][START_REF] Zhao | Direct numerical modeling of time-reversal acoustic subwavelength focusing[END_REF][START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF] with two nonzero surface integrals. The associated concept of the local target model M l1 is depicted in Figure 1d. The hybrid interface S can be freely selected, and the green part corresponds to the PML boundary domain for the local simulation. The required physical hybrid inputs include the gradient ∇q and potential q exactly on the hybrid interface for the acoustic scenarios, as shown in Equation ( 5). The implementation is similar to the VM method but with the full imposing of two surface integrals in Equation ( 6). In the local target model M l1 , the hybrid wavefield can be calculated using the hybrid numerical simulation:

M 1 Q1 + K 1 Q 1 = F T 0 + F q 0 , (11)
where Q 1 is the discrete potential vector, M 1 is the global mass matrix, and K 1 is the global stiffness matrix assigned to the local target model M 1 . F T 0 and F q 0 are the imposing hybrid source terms, which can be constructed using Equations (A2), (A3), and (A4)

with the gradient ∇q and potential q of integration points on a closed arbitrary hybrid interface S, as plotted in the dashed cyan line or red dashed circle in Figure 1d, similar to the red dashed circle in Figure 1a-b. In this study, we selected the GLL points exactly on a rectangular hybrid interface S (the cyan dashed line surrounding the white box depicted in Figure 1d) as the numerical integration points, whose coordinates were determined using only the meshing of the local reference model M l0 .

Through the theoretical proof, for a given closed local domain Ω l , the derivation of Equation ( 20) in [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF] indicates that the two surface integrals on the hybrid interface S in the RP hybrid method will have zero contribution (exact cancellation) outside of the local target domain, because u M i in Equation ( 20) is defined in the

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth entire global domain Ω and not the local domain Ω l . Physically, imposing a traction boundary condition such as Equation ( A3) is equivalent to applying forces on the hybrid interface, whereas imposing a displacement boundary condition such as Equation ( A4) is equivalent to applying moment displacement tensor sources along the hybrid interface [START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF].

The superimposed contribution of these two surface integrals with the ± symbols is +100% inside (full recovery) and zero outside (exact cancellation). [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF] numerically proved it using the FDM in a local reference model M l0 with the spherical shaped hybrid interface, and they indicated that the contributions of the traction and displacement surface integrals are +50% and -50% outside the local domain Ω l and +50%

and +50% inside the local domain Ω l , respectively. Slightly different from their conclusion, in Figure 3, we numerically prove this using the SEM in a local reference model M l0

but with a rectangular shaped hybrid interface. In the global domain, we observe that the two waves in Figures 3a and 3b

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Merits and Limitations

From the above analysis, multiple point sources methods such as the VM and RP methods only require the respective physical hybrid inputs of some integral points (e.g., GLL points around the local mesh) exactly on the hybrid interface S, which is memory- tions mostly rely on degrees 4-8 in each direction, higher degrees are often not adopted, primarily owing to the explicit meshing of mechanical discontinuities and exceedingly small available time steps. Note that in the recent homogenization method in seismology to smoothen the internal mechanical discontinuity [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF], the smooth models used in forward/backward simulations using FWI [START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF], the computational complexity analysis of code-independent features for SEM, and the actual computation time benchmarks all make very high polynomial degrees SEM attractive and competitive [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF].

For the VM hybrid method, the adopted absorbing boundary condition (ABC) is effective in the local target model M l1 with small structural perturbations [START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF]Tong, Chen, et al., 2014), but it is not sufficiently accurate for the scattered tangential incidence waves in complicated local target models in the hybrid calculation [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF][START_REF] Xie | Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML[END_REF]. It cannot adopt the PML absorption because F T 0l can only be imposed on the boundary of the local domain ∂Ω but not freely inside the local target model M l1 in Equation ( 11) owing to the required doubling amplitude of the free reflection from the hybrid interface.

Although the RP method can be naturally combined with the PML absorbing condition owing to the full consideration of the representation theorem, the resultant hybrid wavefield is intrinsically inaccurate in the hybrid interface elements (red elements in Figure 4a1) in which the hybrid inputs are imposed. This is because the wavefields in the elements containing the sources are not accurate for SEM [START_REF] Nissen-Meyer | A twodimensional spectralelement method for computing sphericalearth seismograms i. momenttensor source[END_REF]. In the framework of normal numerical simulations using SEM, recent effective sources can be used to address this problem [START_REF] Capdeville | Homogenization of seismic point and extended sources[END_REF]. In the hybrid simulation framework, an additional layer of elements can be used to bypass its influence [START_REF] Lin | Highfrequency seismic wave modelling of the deep earth based on hybrid methods and spectral-element simulations: a conceptual study[END_REF].

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Direct Discrete Differentiation Method

The second major category, namely, the direct discrete differentiation method, is based on the numerical representation theorem and focuses on the discrete wave equation, primarily consisting the "BY" (abbreviation expression from [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF]), "YM" (abbreviation expression from Masson et al. ( 2014)), and "MYM" (the modified YM) hybrid methods.

BY Hybrid Method

Unlike the multiple point sources method, which is explicitly based on the representation theorem, [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF] and [START_REF] Yoshimura | Domain reduction method for three-dimensional earthquake modeling in localized regions, Part II: verification and applications[END_REF] directly operated on the discrete wave equation to obtain the numerical hybrid inputs using the FEM in the global reference model M g0 . Hereafter, we call it the BY hybrid method, and the associated concept of the local target model M l1 is depicted in Figure 1e. The detailed formulae of the equivalent seismic forces (similar to the numerical hybrid inputs) in the SEM framework are as follows (refer to [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF] and [START_REF] Yoshimura | Domain reduction method for three-dimensional earthquake modeling in localized regions, Part II: verification and applications[END_REF] for more details):

F eff 0 = -K se 0 Q e 0 + K es 0 Q s 0 , (12)
where K is the stiffness matrix, the subscript 0 means that the relevant calculations are performed in the global reference model M g0 , and the superscripts e and s represent GLL points located inside the external domain Ω e and on the hybrid interface S, respectively.

The key property of this numerical hybrid inputs F eff 0 involve only the submatrices K se 0 and K es 0 , which vanish everywhere except in a single layer (the blue part in Figure 1e and the blue elements in Figure 2) in the external domain Ω e adjacent to the hybrid interface S. Thus, the wavefield required to calculate F eff 0 is associated with all the GLL points in the single-layer elements (referred to as the "hybrid elements"), but not only on the hybrid interface S.

The BY hybrid method is an innovative and attractive method in that a teleseismic event can be replaced equivalently by forces loaded on all the GLL points in the hybrid elements in which the numerical hybrid inputs are imposed. Note that in the BY hybrid simulation, after imposing the numerical hybrid inputs F eff 0 , the remaining wavefields going outside owing to the existence of unknown heterogeneities in the local target model M l1 can be naturally eliminated using a suitable absorbing boundary condition (e.g., the PML boundary domain corresponding to green elements in Figure 1e). This numerical expression can obtain numerical hybrid inputs from a 3D background model for 3D local imaging with a higher resolution, namely, box tomography. However, the expression of the hybrid inputs in Equation ( 12) implies storing the global potentials of

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth all the GLL points in the single layer elements, and this requires an excessive disk usage, particularly for 3D numerical simulations [START_REF] Yoshimura | Domain reduction method for three-dimensional earthquake modeling in localized regions, Part II: verification and applications[END_REF], and it might be even worse if people use the very high polynomial degrees SEM, as in [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF].

Note that physical quantities such as the potential/displacement of the GLL points in the hybrid elements can be interpolated from the global mesh to the local mesh, but F eff 0 cannot be precomputed in the global simulation and interpolated into a mesh different from the global mesh, because the internal force term KQ is discontinuous (not smooth)

between connected elements and along the hybrid interface. with fewer storage capacity of numerical hybrid inputs than the outer scheme, can be used to construct the numerical hybrid inputs using

YM Hybrid Method

F M 0 = e W e • f e -W e • (K e 0 • Q e 0 ) + K e 0 • (W e • Q e 0 ) = W • f -W • (K 0 • Q 0 ) + K 0 • (W • Q 0 ). ( 13
)
In contrast, in the global elements crossed by hybrid interface S, in which most GLL nodes belong to the external domain Ω e , the outer scheme can be used to construct the numerical hybrid inputs using

F M 0 = e W e • f e + (I -W e ) • (K e 0 • Q e 0 ) -K e 0 • [(I -W e ) • Q e 0 ] = W • f + (I -W) • (K 0 • Q 0 ) -K 0 • [(I -W) • Q 0 ]. ( 14
)
Where

Accepted Article

This article is protected by copyright. All rights reserved. For the YM method, on the outside of the hybrid interface S, an additional layer of elements is always required to impose the hybrid inputs (e.g., Figure 12b main Ω l , including the hybrid interface S, can be recovered after imposing the hybrid inputs, but only the remaining scattered waves are simulated in the hybrid elements. Note that in the best scenario, as discussed in Section Appendix C, the numerical hybrid inputs of YM method are the same as those of the BY method, but their proposed methods of imposing the hybrid inputs are different. Less memory is required for the YM method because it is only related to the internal force W e (K e 0 •Q e 0 ) and potential W e Q e 0 of the GLL points exactly on the green line in Figure 2; not the Q e 0 of all the GLL points in all the hybrid elements must be stored as suggested for the BY method. However, when the local meshing differs from the corresponding global meshing, spatial interpolation is used to obtain all the displacements of the GLL points required in the hybrid elements during the first global simulation.

MYM Hybrid Method

However, we observe that in the YM method, the additional layer of elements sur-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth the MYM hybrid method, and the concept of local target model M l1 associated with the MYM method is depicted in Figure 1f. Note that the hybrid elements of the MYM method are located inside the local domain, but the hybrid elements of the YM method are outside the local domain selected using different window functions, as indicated by the red and blue elements, respectively, in Figure 2. We can utilize the internal forces (I-W e )•

(K e 0 •Q e 0 ) and (I-W e 0 )•Q e 0 at the GLL points exactly on the hybrid interface S in the hybrid elements to implement the hybrid simulation, according to the outer scheme of the YM method with a new Heaviside window function, depending on the equation ( 14).

It is not an inner scheme compared to the YM method even the hybrid elements are inside the local domain, but more like an outer scheme of the YM method with a new Heaviside window function. Fewer elements can be used to perform hybrid numerical simulations without overlapping the scattered and full wavefields. Note that the red part in The resultant hybrid waveforms are perfectly accurate even when enlarged by the factor of 10 40 when the local meshing is the same as the global meshing surrounded by the green line as shown in Figure 2. The reason for the preciseness is that the MYM hybrid method depends on the numerical hybrid inputs to implement the local hybrid simulation. In the local simulation using the MYM method, all the potentials of the GLL points required by the equation ( 14) are calculated and saved during the global simulation, then a new local equation fully equivalent to the global equation inside the local reference model is compactly reconstructed. A similar numerical approach operating on the matrices of the solvers is shown in [START_REF] Bielak | Domain reduction method for three-dimensional earthquake modeling in localized regions, part i: Theory[END_REF]. We fetch the required local part of the global matrix solver; thus, they should be exactly the same. However, when the local meshing is different from the corresponding global meshing, as in the YM method, we must also use the spatial interpolation to obtain all the potentials of GLL points in the hybrid elements during the first global simulation. The related spatial-interpolation error will affect the results of the hybrid simulations to an extent.

Merits and Limitations

From the above analysis, the BY, YM, and MYM methods have certain advantages:

(i) In the best scenario, only the numerical hybrid inputs of GLL points exactly on the

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth hybrid interface S are required, which is also memory-saving and very suitable for 3D hybrid simulation and is the same as the multiple point sources method. (ii) The PML can be applied naturally to the remaining scattered wave going outside, which results in an accurate hybrid simulation. The MYM method does not require an additional layer of elements outside the local target domain Ω l . However, they are requried by the BY and YM methods to impose the hybrid inputs and with only the scattered waves inside.

Therefore, for the MYM method, the scattered and full wavefield do not co-exist.

However, the direct discrete differentiation method has a disadvantage. 1D global solvers such as the DSM tending to offer physical hybrid inputs without spatial interpolation errors are relatively difficult to apply with the BY, YM, and MYM methods because all the displacement at the GLL nodes in the hybrid elements are required to calculate the internal force term. For 3D global solvers such as SPECFEM3D globe that can offer numerical hybrid inputs, for an accurate implementation of numerical hybrid inputs, the meshing of the local target model M l1 connected with the hybrid interface should be exactly the same as the meshing of the global reference model M g1 . However, when the local meshing is different from the global meshing, we must calculate all the corresponding displacements of all GLL points in the hybrid elements of the local meshing using spatial interpolation and then obtain part of the hybrid inputs by calculating the tensor product KQ inside the hybrid elements. These hybrid methods are all based on the numerical representation theorem; they require a large amount of storage, and they are not suitable for the efficient very high polynomial degree spectral element method.

Because numerous GLL points of the hybrid elements are required to construct the numerical hybrid inputs (approximately N +1 2 times memory as the multiple sources method, N is the polynomial order of the adopted SEM; refer to the discussion section for detailed derivation), and the spatial interpolation introduces some errors, limiting its use to an extent.

New Hybrid Method

As the above analysis shows, the RP method can naturally adopt the PML, and their meshing connected with the hybrid interface of the local and reference model can differ from the global meshing, resulting in the possibility of adopting the more efficient where the hybrid inputs are imposed. The MYM method is capable of natural PML absorption, but the memory requirement is significantly larger than the multiple point sources method particularly for the very high polynomial degrees SEM. For the VM method, the local meshing is independent of the global meshing, but it cannot adopt efficient PML

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth absorption. To obtain the essence of fewer memory requirements, flexible meshing, and natural PML absorption, and discard the dregs, such as only absorption with ABC or inaccurate simulation of the elements with sources, we develop a new hybrid method by combining the multiple point sources and direct discrete differentiation methods.

First, we combine the VM method in Equation ( 10) and MYM method in Equation ( 14).

F T 0 = M 0 Q0 + K 0 • Q 0 F M 0 = Wf -(I -W) • (K 0 • Q 0 ) -K 0 • [(I -W) • Q 0 ]. ( 15
)
After replacing the internal force term K 0 •Q 0 in the second equation with the first equation, we can obtain a new formula for the numerical hybrid inputs:

F M 0 = Wf -(I -W) • (F T 0 -M 0 Q0 ) -K 0 • [(I -W) • Q 0 ]. ( 16
)
We circumvent the requirement for q at all GLL points in the hybrid elements to calculate the internal force term and instead use the difference between the external force term (traction term F T 0 ) and the inertial force term (acceleration term M 0 Q0 ). Therefore, we can utilize only three physical terms on the hybrid interface to obtain hybrid inputs. Using Equation ( 16), to obtain the new combined hybrid inputs F M 0 , we first require the gradient ∇q to construct the traction vector F T 0 , as mentioned in Section 3.5.1.

In addition, we require the potential q and the second temporal derivative of the potential ∂ tt q to build the numerical hybrid inputs during the first global simulation in the reference model M g0 . Note that in practice, we can compute q directly from the saved q field instead of saving them during the global simulation. These three physical quantities are located at the integration points (e.g., GLL) exactly on the hybrid interface S. To implement the new hybrid simulation, we can directly impose the hybrid inputs into the integration points on the hybrid interface regardless of whether the internal structure is changed. We can conclude that the proposed hybrid method has the following advantages:

• The three physical quantities used to construct the hybrid inputs can be obtained flexibly using the global solver (e.g., DSM, SPECFEM3D globe, or AxiSEM3D), and the meshing of the local target model is completely independent of the meshing of the global reference model.

• Only the associated three physical quantities of the points exactly located on the 2D hybrid interface are required, which is highly accurate and memory efficient

-19-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth for the 3D hybrid numerical simulation and can increase the computational efficiency of the box tomography.

• The efficient PML can be adopted naturally without any elements overlapping between the local and PML domains, which can further increase the imaging accuracy of the box tomography.

The concept of the local target model M l1 associated with the new hybrid method is depicted in Figure 1f. Figures 4d1 and4d2 show the hybrid wavefield and waveforms of the new hybrid method, and we can observe that the error of the new method is almost the same as the VM method in Figure 4b2 for the same local reference model M g0 . However, as discussed in Section 3.5.1 and the numerical simulations benchmark in Section 4. 1.3, if the local target model is strongly heterogeneous, the VM method is not effective because it cannot be combined with the PML absorbing condition. The RP hybrid method would be a better alternative compared with the proposed method when the traditional SEM with NGLL = 5 is used for the box setting with the source outside and receiver inside (SORI). However, if the very high polynomial degree SEM is used in the hybrid simulation, and the acceleration q is calculated using the stored potential q during the hybrid simulation, the new method will have a more efficient computation cost of the local simulation and fewer memory requirements of the hybrid inputs. Note that combining the new and RP hybrid methods is very useful for the hybrid simulations in the box settings with the source inside and receiver outside (SIRO) and the source outside and receiver outside (SORO). Note that a special hybrid method for separately solving the PML domain and full wavefield domain was implemented in [START_REF] Xie | Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML[END_REF]. Their key concept of their hybrid method was to separate the entire domain into the diffracted waveform domain and full waveform domain, and to solve them separately in two different ordinary differential equations, while our proposed method solved them in one ordinary differential equation same as the traditional solver using the spectral element method.

Note that the RP hybrid method theoretically results in zero contributions on the hybrid interface S when the hybrid simulation is performed in the local reference model M l0 . We must consider the known saved potentials and the scattered potentials obtained by the hybrid simulation on the hybrid interface when calculating the hybrid waveforms of the stations located within the elements containing the hybrid interface. It is the same for the MYM and "new" hybrid methods, but it is not required for the VM, BY, and YM hybrid methods.

In summary, in this section, we analyze in detail the two main hybrid categories and five subcategories (RP, VM, BY, YM, and MYM) of hybrid methods, and we propose a new hybrid method. The new hybrid method inherits the advantages of the VM

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth and MYM methods and eliminates their disadvantages. In Table 1, we compare these six hybrid methods in four criteria under 3D hybrid simulations as in Sections 4.3 and 4.4: i) physical or numerical representation theorems; ii) stored quantities in the best scenario or when the local and global meshings are different; iii) the required storage size in a 3D hybrid simulation in the best scenario or when the local and global meshings are different; iv) whether the absorbing boundary condition can be used. In the next section, we discuss a series of numerical experiments to benchmark and verify these hybrid methods.

Numerical Experiments

To numerically validate the proposed hybrid method and benchmark its numerical accuracy with other existing hybrid methods, we set up a series of 2D and 3D numerical experiments in homogeneous and heterogeneous acoustic models. For the homogeneous models, we utilized V P 0 = 3750 m/s for the P-wave velocity and ρ 0 = 2000 kg/m 3 for the density. For the heterogeneous scenario, we kept the density constant, and κ(x)

was in the form of a Gaussian or monochromatic oscillatory function. The detailed definition of heterogeneity was defined by a spatial function f (x) as follows:

κ(x) = κ 0 f (x), (17)
where

κ 0 = ρ 0 V 2 P 0 .
For the descriptions of the model size and the time duration in the acoustic scenario, the maximum frequency and P-wave velocity defined the minimum wavelength of the propagation wavefield as

λ min = V P f max . ( 18
)
Subsequently, we measured the spatial scale as a function of λ min and the temporal scale as a function of

t min = 1 f max . (19)
The source was a point force located in the middle of the free surface, and its time wavelet was a Ricker function (the second derivative of a Gaussian) with a central frequency of f 0 = 2 Hz for 2D scenarios and f 0 = 0.4 Hz for 3D scenarios. Thus, λ min was approximately 625 and 3125 m in 2D and 3D scenarios, respectively. Thirteen receivers were evenly located in the local models, as shown in Figures 2 and10a for the 2D and 3D scenarios, respectively. In the global numerical simulations, the boundary conditions were free normal stress conditions all around the domains to fully exclude imperfect absorption, and the selected domains were sufficiently wide to ensure that no reflected waves

Accepted Article

This article is protected by copyright. All rights reserved.

N GLL points in each direction, nax0 = (nex -2) × N -1, nay0 = (ney -2) × N -1, naz0 = (nez -2) × N -1; nax1 = nex × N + 1, nay1 = ney × N + 1, naz1 = nez × N + 1; nax2 = nex × N -1, nay2 = ney × N -1, naz2 = nez × N -1; nax3 = (nex + 2) × N + 1, nay3 = (ney+2)×N +1, naz3 = (nez+2)×N +1; nax4 = (nex+2)×N -1, nay4 = (ney+2)×N -1, naz4 = (nez + 2) × N -1. Where N = N GLL -1. The number of GLL points in a closed 2D hybrid interface is V1 = (nax1 × nay1 × naz1 -nax2 × nay2 × naz2).
When the local and global meshings are different, the number of GLL points used to construct the numerical hybrid inputs using the MYM method is V2 = (nax1 × nay1 × naz1 -nax0 × nay0 × naz0); the number of N GLL points used to construct the numerical hybrid inputs using the YM method is

V3 = (nax3 × nay3 × naz3 -nax2 × nay2 × naz2)
, and the number of GLL points used to construct the physical hybrid inputs using the RP method is V4 = (nax3 ×nay3 ×naz3 -nax4 ×nay4 ×naz4), because one layer of elements inside the hybrid interface are needed to ensure simulated accuracy.

nt is the number of time steps in the simulation without temporal interpolation. For example, the storage size of the VM is 2 × 3V1 × nt, where 2 represents the items including traction and velocity, and 3 indicates three directions in 3D. For the absorbing condition, ABC is a type of absorbing boundary-type condition, e.g., [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF]. Note that for the new method, q is calculated directly from the saved q field to decrease the amount of storage.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth from boundaries affected the results. The simulated duration and temporal step ∆t were 90 t min (15 s), 0.002 s for 2D scenarios, and 24 t min (20 s), 0.01 s for 3D scenarios.

One of the key parameters when using the SEM is G (the number of grid points per minimum wavelength). Classic SEM applications mostly rely on degrees 4-8 spectral elements in each tensorial direction, and the G for a high degree 20 is approximately half of that necessary for a degree 4 [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. The commonly admitted G to obtain sufficient accuracy in a constant-velocity medium is approximately G ≈ 5 for one element with degree 8 (NGLL = 9, number of GLL points per element per direction) [START_REF] Basabe | Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF][START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. A low G can be important in the FWI context owing to the extensive operations on forward and adjoint wavefields [START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion[END_REF][START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF]. In our study, the meshing based on NGLL = 8 was adopted for the first global simulations in the global reference model M g0 . A different meshing with NGLL = 20 was used for the VM, RP, and new hybrid simulations in the local reference models M l0 and target models M l1 on account of the efficiency of very high polynomial degree SEM [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. This article is protected by copyright. All rights reserved.

2D Homogeneous Global Reference Model

manuscript submitted to JGR: Solid Earth imposed was very large (≈ 10%) owing to the intrinsic inaccuracy of the source elements in the SEM (Figures 4a1 and4a2). The relative errors of the VM and new methods were almost the same (≈ 0.01%) (Figures 4b1,4b2,4d1,and 4d2), which originated from the boundary reflection of the spatial-dispersion error difference of the global and local meshing. Through the hybrid numerical simulations using the RP, VM, MYM, and new methods in the same 2D global and local homogeneous reference model, we initially demonstrated that the accuracy of the new method is the same as that of the VM method, but as discussed in Section 3.5.1, if the local target model is strongly heterogeneous, the VM method is not effective because it cannot be combined with the PML boundary condition. Here, the MYM method is more accurate.

Hybrid Simulation for a Local Weak Gaussian Heterogeneous Model

To validate the applicability of the new hybrid method in the 2D weak heterogeneous model, the combined hybrid inputs obtained in the 2D homogeneous reference model were imposed into the local target model M l1 with a Gaussian heterogeneity inside. The adopted 2D Gaussian spatial function f (x) was

f (x) = 1 + ae -x 2 2σ 2 , ( 20
)
where a is the amplitude value that controls the range of κ(x), and σ is the spatial scale parameter. Here, we used a = -0.

Hybrid Simulation for a Local Strong Gaussian Heterogeneous Model

To validate the applicability of the new hybrid method in the 2D strongly heterogeneous model, the hybrid inputs that were obtained in the 2D homogeneous reference model were imposed into the local target model M l1 with a strong Gaussian heterogeneity inside. Here, we used a = -0.8 and σ = 2 λ min to construct a global Gaussian model.

The parameters of the meshing were the same as those in the weak scenario. Figure 6a shows the hybrid wavefield using the VM hybrid method. The hybrid inputs were calculated in the 2D homogeneous reference model. It was clear that the wavefront was severely deformed after passing through the low velocity. The relative waveform error of using the VM method was approximately 7.5% (Figure 6b), and the relative waveform error of using the new method was approximately 0.89% (Figure 6c) but with more calculation time owing to the use of the PML absorbing condition.

To further verify the applied scope of the proposed hybrid method, we further tested it in 2D heterogeneous, 3D homogeneous, and heterogeneous reference models M g0 .

2D Global Heterogeneous Reference Model

For the 2D heterogeneous scenario, the size of the 2D global heterogeneous reference model M g0 was the same as that of the 2D homogeneous scenario. The adopted 2D cosine spatial function f (x) was defined as follows:

f (x) = 1 + a cos( 2π λ h k a • x) + cos( 2π λ h k b • x) k a = cos(t 1 ), sin(t 1 ) -25-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

k b = cos(t 2 ), sin(t 2 ) ( 21
)
where a is the amplitude value, and t 1 , t 2 are the included angles between the two heterogeneities and the x-axis. In this scenario, we used a = 1 20 and t 1 = π 4 , t 2 = 3π 4 (in radians) related to two orthogonal directions (Figure 7a); λ h = λ min corresponds to a typical geological model [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF].

Figure 7b shows the global domain Ω g and the wavefield snapshot at 67.2 t min associated with 240 × 120 elements meshed with NGLL = 8 (G = 12, 1.5 times as the homogeneous scenario based on the empirical value in [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. The local domain Ω l and wavefield snapshot at the same time step associated with 30 × 15 elements meshing with NGLL = 20 (G = 7.5), for the scenario λ h = λ min , are depicted in Figure 7c.

The agreement between the hybrid and reference waveforms from the 13 receivers is shown in Figure 7d, demonstrating the effectiveness of the new hybrid method in the 2D heterogeneous global reference model M g0 .

We know that in actual hybrid applications, there are three types of source and receiver settings: the first one is the box with the SORI, the second one is the reverse with the SIRO, and the third is a box with the SORO. All the above-mentioned 2D hybrid simulations corresponded to the first SORI setting. In the following, we discuss how to perform the other two types of hybrid simulations by adopting the suitable hybrid method(s). ii) perform a global numerical simulation using a single force at the location of the receiver, and record the Green's function and its gradient of the GLL points on the hybrid interface; iii) apply a quadrature rule to discretize the two surface integrals in the representation theorem (Equation 5), perform the convolution at each integral point (e.g., GLL points), and sum the contributions from all the integral points as follows.

q(x, t) = ∞ -∞ dτ ( S f - S ) 1 ρ ∇ q (x , τ ) G (x, t -τ ; x , 0) • ndS(x ) - ∞ -∞ dτ ( S f - S ) 1 ρ q (x , τ ) ∇ G (x, t -τ ; x , 0) • ndS(x ) = 0 -∆t e Ns p=1 α e,p conv( 1 ρ p n j ∂ j q(x p , t), G(x p , t)) -conv( 1 ρ p q(x p , t), n j ∂ j G(x p , t)) . (22)

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Where the first 0 means the zero contribution from the free surface S f (∇ q and ∇ G are equal to zero), the first summation is assembled over all the surface elements e of the hybrid interface S, and the second one is over all the quadrature points x e p (p = 1, • • • , N s ), N s is the total number of quadrature points on one surface element with the related quadrature weight α p , and conv(f (t), g(t)) is the convolution between two time series f (t) and g(t).

The ordering of steps i) and ii) can be changed. It is worth noting that, in this SIRO setting, the stations are always located far from the hybrid interface, and we can explicitly adopt the representation theorem by using the traction and potential terms obtained from the local simulation from the source side and using the Green's function and its gradient obtained from the global simulation from the receiver side. We propose not to adopt the new hybrid method in this SIRO setting.

Figure 8a shows an example of the 2D hybrid wavefield related to the calculation of Green's function at time step 72 t min from the receiver side, corresponding to step ii).

Note that there are some stripes in Figures 8a and9a (setting SORO); these are because the adopted spatial meshing could not accurately simulate the wavefield with the delta source time function for the calculation of the Green's function. Please note that high frequency stripes abovementioned can be implicitly filtered out by the convolution between the Green's function with unlimited frequency band and the forward field with limited frequency band obtained from the banded source. In Figure 8b, the waveform benchmarks are listed. The waveform was explicitly calculated using the convolution of the representation theorem, and it had a very good consistency with the waveform calculated in the global reference model.

Hybrid Simulation with a Source Outside and a Receiver Outside

For a remote local target domain in the SORO setting, the workflow has the following four steps: i) perform a numerical simulation in the global reference model from the remote source side, as in the first step in the SORI setting, and record the traction, potential, and acceleration terms of the GLL points on the hybrid interface required by the new hybrid method; ii) perform a numerical simulation in the global reference model from the remote receiver side using a single force, and record the Green's function and its gradient of the GLL points on the hybrid interface; iii) perform a hybrid numerical simulation using the new hybrid method in the local target model by imposing the hybrid inputs recorded in step i), and record the residual traction and residual potential of the GLL points on the hybrid interface; iv) apply a quadrature rule to discretize the two surface integrals in the representation theorem (Equation 5), perform the convolu-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth tion at each integral point (GLL points), and sum the contributions from all the integral points. Note that the convolutions are explicitly calculated between the residual traction and potential from step iii) and Green's function and its gradient from step (ii). We propose adopting the new hybrid method and representation theorem together in the SORO setting.

Figure 9a shows the 2D hybrid wavefield related to the calculation of Green's function in step ii) at time step 72 t min for the SORO setting. Because the distance between source and receiver was larger, the duration was 17 s (102 t min ) in this scenario. In

3D Global Heterogeneous Reference Model

Finally, we performed a hybrid numerical simulation using the 3D cosine heterogeneous model shown in Figure 11a. The size of the global heterogeneous model M g0 was the same as in the 3D homogeneous scenario. The 3D cosine spatial function f (x) was

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth adopted as follows:

f (x) = 1 + a cos( 2π λ h k a • x) + cos( 2π λ h k b • x) + cos( 2π λ h k c • x) k a =
where a is the amplitude value, and (t 11 , t 12 ), (t 21 , t 22 ), (t 31 , t 32 ) are the included angles between the three heterogeneities along the x-and y-axes, respectively. Here, we used a = 1 30 and (t 11 = 0, t 12 = 0), (t 21 = 0, t 22 = π 2 ), (t 31 = π 2 , t 32 = π 2 ) corresponding to three orthogonal directions shown in Figure 11a, and λ h = λ min was adopted, corresponding to a typical geological model [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF]. The global domain

Discussion

Spatial and Temporal Interpolations

In the 3D numerical simulations of wave propagation, dividing the minimum sampling points per minimum wavelength G by 2 implies a storage requirement of the hybrid inputs divided by 4 because of the 2D planar distribution, which is significant for decreasing the memory requirements. For example, for the spatial interpolation in our 3D homogenous scenario, G = 8 in the local simulation resulted in the storage of the hybrid inputs being 1.6 2 times that of G = 5 in the local simulation (the 3D homogeneous scenario). Moreover, the computation time with degrees 12-24 can be up to twice as fast as the classic degree 4 [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF]. Considering a 3D hybrid local simulation, the corresponding hybrid inputs are a 2D planar distribution; thus, the storage demand of the hybrid inputs of degree 20 (G = 2.5) is a quarter of that of degree 4 (G = 5).

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth hybrid elements, and finally subtracts the calculated internal force from the internal force in the target model based on Equation ( 14). The YM and BY methods use opposite steps based on Equation ( 13). The different hybrid methods were designed to perform the hybrid simulation; however, the numerical complexity of the SEM solvers is primarily dominated by the calculation of the internal forces (the product between the elemental stiffness matrix and displacement vector) [START_REF] Deville | High-order methods for incompressible fluid flow[END_REF], and the imposing of hybrid inputs using different hybrid methods into the local SEM solver does not change the main computation part, because the number of hybrid elements is significantly lower than the total number of elements in the local domain.

Hybrid Inputs Obtained from Programs SPECFEM3D globe and AxiSEM3D for Box Tomography

Box tomography has been used to image the upper-mantle shear velocity structure beneath the North American continent down to a 40 s period [START_REF] Clouzet | Box Tomography: first application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[END_REF]. Their hybrid inputs are calculated by the program SPECFEM3D globe [START_REF] Komatitsch | The spectral-element method, beowulf computing, and global seismology[END_REF] for the seismic events outside the box based on the YM hybrid method. There is nothing that prevents one from using the new hybrid method to obtain the hybrid inputs in the framework of SPECFEM3D globe. Alternatively, with two to three orders of magnitude faster, the program AxiSEM3D [START_REF] Leng | Axisem3D: broadband seismic wavefields in 3-D global earth models with[END_REF]) is an excellent and efficient candidate solver for calculating the hybrid inputs in an existing 3D global reference model.

Subsequently, hybrid simulations can be implemented with a local solver (e.g., RegSEM [START_REF] Cupillard | Regsem: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale[END_REF]). It is also worth emphasizing that the global solver, which is used to provide the hybrid inputs, only requires to calculate the physical quantities displacement, acceleration, and strain (for the elastic scenario) of the integration (e.g., GLL) points located exactly on the hybrid interface. The local solvers used for the hybrid numerical simulations only require to read in these physical hybrid inputs as equivalent virtual sources. There are no other connections between these two programs; therefore, the proposed hybrid method is highly accurate and memory efficient for implementing hybrid simulations in actual applications. Note that the proposed hybrid method, similar to [START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF] and [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF], cannot recover second-or higher-order scattered waves, but first-order scattered waves are always stronger and more important for imaging.

Conclusion

By analyzing and combining the respective advantages of the existing two main categories and five subcategories of hybrid methods, a new hybrid method is proposed. Only

  Figure 2 shows the 2D wavefield in a global reference model M g0 ; note that the snapshots at five different time steps are superposed to simplify the diagram. Figures 4b1 and

  Figure 3a-b can recover the hybrid wavefield, as shown in Figure 3c. Therefore, the RP hybrid method can obtain the full recovery inside the local domain and with a possible external domain, has a natural ability to combine the PML.

  saving and suitable for the practical 3D hybrid simulations. The meshing of the local target model M l1 is completely independent of the meshing of the global reference model M g0 . The associated local hybrid simulation is flexible, resulting in the possibility of using a more efficient SEM with very high polynomial degrees (e.g., degrees 12 to 24)[START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF] to reduce memory requirements, speed up the forward simulation, and use the upscale non-periodic homogenization[START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF][START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF] to improve the ability of FWI. Although classic SEM applica-

  [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF] further developed the BY method and proposed a more compact numerical hybrid method, which can accurately perform the hybrid simulation in the local target model with arbitrary shapes. Hereafter, we call it the YM method, and the associated concept of local target model M l1 is depicted in Figure 1e. Note that only the best scenario (the same meshing used in local and global domain, described in Section Appendix C) is plotted here to illustrate the arbitrary hybrid interface; please refer to Masson et al. (2014) for more details. The YM method can construct the mirror excitation F M 0 (the same as the numerical hybrid inputs) with two different schemes, which relies on the spatial window function and discrete wave equation. For a local target model M l1 with arbitrary shapes, in the global elements crossed by the hybrid interface S, in which most GLL nodes belong to the closed local target domain Ω l , the inner scheme,

e

  denotes the assembly of all element hybrid inputs. Q e 0 , K e 0 and Q 0 , K 0 are the element and assembled potential vectors and stiffness matrices, respectively, assigned to the global model M g0 associated with the global domain Ω g . The element and assembled diagonal matrices W e and W are discrete window functions used to fetch the local wavefield in the local domain from the global wavefield. For example, W e •Q e 0 and -15-

  manuscript submitted to JGR: Solid Earth W•Q 0 are related to the element and assembled potential vectors, respectively, in the local domain Ω l . The diagonal entries of W e = 0 represent the GLL points in the global elements crossed by the hybrid interface S, and they simultaneously belong to the local target domain Ω l . The diagonal entries of W e = 0 represent all the other GLL points. Both the inner and outer schemes must record the local internal force W e • (K e 0 • Q e 0 ) or the complement (I -W e ) • (K e 0 • Q e 0 ) and the potential W e • Q e 0 or the complement (I-W e )•Q e 0 in the same hybrid elements crossed by the hybrid interface S during the first global simulation. The only elements that contribute to the numerical hybrid inputs F M 0 are the hybrid elements, for which the value of the window function W e evaluated at the GLL nodes is not constant. When the hybrid inputs are imposed into the blue part (the hybrid elements in SEM) shown in Figure 1e, the inner scheme adds the saved displacement first before calculating the internal force, and then subtracts the saved internal force, while the imposition of the outer scheme subtracts the displacement first and then adds the internal force (refer to Masson et al. (2014) for the detailed derivation).

  [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF] and Figure5in[START_REF] Clouzet | Box Tomography: first application to the imaging of upper-mantle shear velocity and radial anisotropy structure beneath the North American continent[END_REF]). The full wavefield in the local do-

  rounding the local target model (the blue region in Figure1eand blue elements in Figure 2) is not necessary. Analogously, a new Heaviside window function can be adopted such that the diagonal values W e = 1 are only in the local domain Ω l , but W e = 0 are both on the hybrid interface S and in the external domain Ω e . Hereafter, we call it

Figure

  Figure 1f and red elements in Figure 4c1 are the hybrid elements that are selected using the new window function W e to locate inside the local target domain, and the green part corresponds to the PML domain. In the best scenario, the associated numerical hybrid inputs are only imposed on the GLL points exactly on the cyan dashed line in Figure 1f and on the green line in Figure 4c1 between the PML and hybrid elements.

Figures

  Figures 4c1 and 4c2 depict the hybrid wavefield and waveforms of the MYM method.

For

  2D homogeneous scenarios, the size of the global reference model M g0 was 160 × 80 λ 2 min , which was assigned to a rectangular global domain Ω g . The domain Ω g and wavefields overlapping at five time steps associated with the 160 × 80 structural spectral elements with NGLL = 8 (G = 8) are plotted in Figure 2. 4.1.1 Hybrid Simulation for a Local Reference Homogeneous Model The size of the adopted local reference model M l0 was 80 × 40 λ 2 min , which was assigned to the rectangular local domain Ω l surrounded by the closed green line in Figure 2 (hybrid interface S), and it was located at the center of the global reference model. The related wavefield snapshots at the same time steps as the global simulation associated with 20 × 10 elements with NGLL = 20 (G = 5) are shown in Figures 4a1, b1, and d1, which correspond to the RP, VM, and new hybrid methods, respectively. They all had the characteristics of flexible meshing, namely, the number of elements or/and the internal degree of the local meshing for the hybrid simulation differed from the related global meshing of the local domain Ω l . In comparison, the MYM method adopted local meshing with 80 × 40 elements with NGLL = 8 (G = 8), which was identical to the global meshing (Figure 4c1). A comparison of wavefields and waveforms indicated that the MYM method had no error relative to the global simulation shown in Figures 4c1 and4c2. The relative error of the RP method in the elements in which hybrid inputs were

  2 and σ = 2 λ min to construct a global Gaussian model, and the Gaussian heterogeneity existed only in the local domain Ω l . The external model M e was the same as the homogeneous reference model M g0 (Figure 5a).

Figure

  Figure5bshows the hybrid wavefield by imposing the hybrid inputs calculated in the 2D homogeneous reference model into the global Gaussian (target) model M g1 . The hybrid inputs can be considered new equivalent sources to replace the remote source to obtain the same wavefield in the local domain. Note that the meshing used in Figure5bwas constructed using 64 × 32 elements with NGLL = 20 (G = 5), different from the global meshing (160 × 80 elements with NGLL = 8 (G = 8)), which was used for calculating the global waveforms in the Figure5d. In this scenario, the remaining scattered waves going outside were observed owing to the presence of the local Gaussian heterogeneous anomaly. Because these waves were not part of the original reference wavefield,

4. 2 . 1

 21 Hybrid Simulation with a Source Inside and a Receiver Outside For a local hybrid simulation in the SIRO setting, the workflow has three steps: i) perform a local numerical simulation by loading the internal source in the local target model, and record the physical quantities gradient and potential terms (physical quantities needed by the representation theorem) of the GLL points on the hybrid interface;

  Figure 9b, the waveforms benchmark are listed. The black line is the waveform difference between the global waveforms respectively calculated in the global weak Gaussian model in Figure 5 and the reference model in Figure 2, and it had a very good consistency with convoluted waveform. Note that the residual traction and potential were nonzero owing to the local Gaussian heterogeneity. It is also worth noting that in the SORO setting when the local model contained the nonzero local perturbation, the combined contribution of the two surface integrals to the waveform of the remote receiver were nonzero, and it is only the scattered waveform; and when the local target model was the same as the local reference model M l0 without any local perturbation, the combined contribution of the two surface integrals to the waveform of the remote receiver were zero, not the full waveform.4.3 3D Global Homogeneous Reference ModelFor the 3D homogeneous scenario, the size of the 3D global homogeneous reference model M g0 was 32 × 32 × 16 λ 3 min . Figure10ashows the global domain Ω g and the wavefield snapshot at time step 15.6 t min associated with 32 × 32 × 16 elements meshed with NGLL = 8 (G = 8). The local domain Ω l for the hybrid simulation was also a parallelogram located at the center of the global reference domain with a half-size of 16 × 16 × 8 λ 3 min . The hybrid wavefields at the same time step associated with 4 × 4 × 2 elements with NGLL = 20 (G = 5) are depicted in Figure 10b. The agreement between the hybrid and reference waveforms from the 13 receivers is depicted in Figure 10c, initially presenting the effectiveness of the new hybrid method in the 3D global homogeneous reference model M g0 .

Ω

  g and the wavefield snapshot at time 15.6 t min associated with 48 × 48 × 24 elements meshed with NGLL = 8 (G = 12, 1.5 times of the 3D homogeneous scenario based on the empirical value in Lyu et al. (2020)) are depicted in Figure 11b. The local reference model M l0 was located at the center of the global domain Ω g with a half-size of 16 × 16 × 8 λ 3 min . The hybrid wavefield at the same time step and the local domain Ω l associated with 6 × 6 × 3 elements with NGLL = 20 (G = 7.5, 1.5 times the 3D homogeneous scenario) are depicted in Figure 11c. The agreement between the hybrid and reference waveforms shown in Figure 11d demonstrated the effectiveness of the new hybrid method in the 3D global heterogeneous reference model, indicating a significant application prospect in probing the key 3D structures in the deep earth using box tomography.

Figure 1 .

 1 Figure 1. Nomenclatures of the hybrid simulation. (a) The global reference model Mg0 contains the gray external model Me and the white local reference model M l0 , which are assigned to the global domain Ωg, external domain Ωe, and local domain Ω l , respectively. The red and cyan dashed lines represent hybrid interfaces S. (b) The global target model Mg1 contains the gray external model Me and the white local target model M l1 , which are also assigned to Ωg, Ωe, and Ω l . The gray cartesian circle represents local target heterogeneity. (c) Local target model M l1 and its hybrid interface S (cyan dashed line) of the VM hybrid method. (d) Local target model M l1 and its hybrid interfaces S (red and cyan dashed lines) of the RP hybrid method, and the green part represents the PML domain. (e, f) The local target model M l1 of the BY, YM, MYM, and new hybrid methods with the green PML part the hybrid elements in SEM (blue part for BY and YM methods, red part for MYM and new methods), and their hybrid interfaces S (cyan dashed lines).

Table 1 .

 1 Benchmarks of the VM, RP, BY, YM, MYM, and new hybrid methods. For a given 3D local target model buried in the deep earth with nex, ney, and nez elements with

			manuscript submitted to JGR: Solid Earth
		Representation	Stored quantities	Required storage size	Absorbing
		theorem	same/different meshing	same/different meshing	condition
	VM	physical	traction+velocity	2×3V 1 × nt	ABC type
			same	2×3V 1 × nt
	RP	physical	traction+potential	2×3V 4 × nt	any
			same	2×3V 4 × nt
	BY	numerical	potential	3V 3 × nt	any
			same	3V 3 × nt
	YM	numerical	internal force+potential	2×3V 1 × nt	any
			only potential	3V 3 × nt
	MYM	numerical	internal force+potential	2×3V 1 × nt	any
			only potential	3V 2 × nt
	NEW	combination	traction+potential (acceleration)	2×3V 1 × nt	any
			same	2×3V 1 × nt

  cos(t 11 )cos(t 12 ), cos(t 11 )sin(t 12 ), sin(t 11 ) k b = cos(t 21 )cos(t 22 ), cos(t 21 )sin(t 22 ), sin(t 21 ) k c = cos(t 31 )cos(t 32 ), cos(t 31 )sin(t 32 ), sin(t 31 ) ,

they were not cancelled out when they crossed the hybrid interface S, and some residual wavefield leaked out of the local target model M l1 (Figure 5b). We could shrink our model and utilize the PML to absorb the scattered waves (Figure 5c). In the hybrid simulations shown in Figures 5b and5c, after the wavefront across the Gaussian heterogeneity with lower velocity than the reference, it propagated more slowly than the reference wavefront in the 2D homogeneous model. Additionally, it could be observed from the -24-

Accepted Article

Acknowledgments First, the authors would like to acknowledge Barbara Romanowicz and Sevan Adourian from the University of California, Berkeley, for their comments on the detailed discussions on box tomography. The authors are grateful to all the participants in the HIWAI project for many discussions on computational seismology. The authors also thank the working environment support during the writing from the Beijing Baijiatuan Earth Science National Observation and Research Station, Beijing, China, 100095. We appreciate the editors Michael Bostock and Sebastien Chevrot, and the reviewer Qinya Liu and the anonymous reviewer for their useful comments that improved the manuscript. This study was supported by the National Nature Science Foundation of China (grant nos. 41625016, 41888101, and 42004045), the China Scholarship Council (File No. 201804910289), and the French National Research Agency (ANR-16-CE31-0022-01, 2016-2021). We note that there are no data-sharing issues because all of the numerical parameter information is provided in the figures produced by solving the global and local hybrid numerical simulations in this study. Data are available at the following site: https://doi.org/ 10.6084/m9.figshare.14559666.v1.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth A low G can also be very important in the FWI context because of the kernel operation on the partial, compressed, or full storage of the forward and adjoint wavefields [START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion[END_REF][START_REF] Boehm | Wavefield compression for adjoint methods in full-waveform inversion[END_REF][START_REF] Fichtner | Full seismic waveform tomography for uppermantle structure in the australasian region using adjoint methods[END_REF][START_REF] Lyu | Intrinsic non-uniqueness of the acoustic full waveform inverse problem[END_REF]. For the proposed hybrid method, during the first global simulation, the physical quantities ∇q, q, and ∂ tt q are calculated and recorded to construct the combined hybrid inputs for the hybrid simulation. The flexible meshing in the local target model can differ from the meshing of the global model, resulting in the possibility of using a low G of the hybrid inputs and the associated hybrid simulation with a very high polynomial degree SEM. When the local and global meshings are different, spatial Lagrange interpolation can be adopted to obtain the physical quantities of GLL points exactly on the hybrid interface during the first global simulation. If the global or local meshing is not sufficiently accurate, different spatial dispersion errors will be introduced into the local simulation, resulting in inaccurate waveforms. Note that physical quantities such as the potential/displacement of the GLL points for constructing the hybrid inputs can be interpolated from the global to the local mesh. However, the internal force term KQ cannot be pre-computed in the global simulation and interpolated on a different mesh from the global mesh, because the internal force term is discontinuous between hybrid elements. Note that [START_REF] Monteiller | On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF] also describes that the spatial and temporal interpolation schemes are key to increasing the efficiency of the algorithms in the global hybrid simulation.

In addition to the spatial interpolation using different local meshing with a low G, we can further adopt different temporal interpolations. The hybrid inputs of the proposed hybrid method are entirely based on the physical quantities ∇q, q, and ∂ tt q of the integration points (for example, here we used GLL points), which are localized exactly on the hybrid interface S of the local target model M g1 . During the first global simulation, we should be able to record them at the Nyquist sampling time steps, and then the physical quantities can be recovered (interpolated) according to the time step dynamically determined using the local structures after imposing the recorded physical quantities and before starting the hybrid simulation. Thus, combining the spatial and temporal interpolations should significantly reduce the memory required for hybrid numerical simulations by several tens of orders of magnitude. If different global and local time steps are adopted, different temporal dispersion errors will be introduced into the local simulation, resulting in inaccurate waveforms to an extent. Note that the spatial and temporal dispersion errors have been proven to be irrelevant [START_REF] Koene | Eliminating time dispersion from seismic wave modeling[END_REF][START_REF] Lyu | Removing the Courant-Friedrichs-Lewy stability criterion of the explicit time-domain very high degree spectral-element method with eigenvalue perturbation[END_REF], which provides a theoretical basis for eliminating the temporal dispersion errors in the global and hybrid simulations.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Numerical cost

To use the proposed method, we must know three incident fields (the gradient, potential, and acceleration), with the advantage of being able to use the PML in the scenario of very heterogeneous media without much artificial error owing to incomplete absorption. From the benchmark of different hybrid methods in Table 1, the storage between the new hybrid method and VM or RP hybrid methods is the same whether the local and global meshings are the same, because the q can be obtained from the stored q field. Generally, the local and global meshings are always different because the local and global solvers are different, and the corresponding storage ratio between the new hybrid method and BY or YM hybrid methods is expressed as follows:

The ratio between the new hybrid way and MYM hybrid methods is expressed as follows:

because V 2 ≈< (N +1)×V 1 . Thus, for the very high polynomial degrees SEM such as N = 20 in the 3D scenarios, the proposed hybrid method requires significantly fewer memory requirements than the BY, YM, and MYM hybrid methods.

To build a clear work flow, we always must i) determine the coordinates of the hybrid GLL points in the local domain required using different hybrid methods; ii) calculate and record the physical quantities of the recorded hybrid GLL points during the first global forward simulation, such as the potential/displacement for YM method; iii) impose the hybrid inputs into the local model to perform the local hybrid simulation at each time step. According to the expressions of different hybrid inputs of the six listed hybrid methods (Table 1), the methods of imposing the hybrid inputs are different. At each time step, the VM method explicitly implements the traction surface integral with the integrated GLL points on the hybrid interface based on Equation (A3). The RP method requires an additional term by further explicitly implementing the second potential/displacement surface integral based on Equation (A4). The new method first implements the traction surface integral as the VM method, then directly add the potential/displacement before calculating the internal force terms, and finally directly subtracts the recorded reference acceleration from the acceleration in the target local model based on Equation ( 16). The BY, YM, and MYM hybrid methods are performed almost the same. The MYM hybrid first directly adds the potential/displacement to the hybrid GLL points on the hybrid interface before calculating the internal force, then calculates the internal force using the elemental stiffness matrix and the potential/displacement of all the GLL points in the

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth three physical quantities, including ∇q, q, and ∂ tt q, must be obtained during the first global simulation to construct the hybrid inputs, and the meshing of the local target model is not limited by global meshing. Furthermore, only the associated physical quantities of the integration points that are located exactly on the hybrid interface are required, which is highly accurate and memory efficient for the 3D local hybrid simulation. The required memories of hybrid inputs can be further reduced fourfold if the very high polynomial degree spectral element method is used for the local hybrid simulation. Finally, the more efficient absorbing boundary condition PML can be adopted naturally without any elements overlapping between the local and PML domains. All these characters can increase the imaging resolution and accuracy of box tomography. After theoretical analysis and numerical experiments, the consistencies of the wavefield and waveform between the global and hybrid numerical simulations in 2D/3D homogeneous and heterogeneous models validate our proposed hybrid method, indicating a significant application prospect in probing the key 3D structures in the deep Earth using box tomography.

It should also be noted that this study focused only on the acoustic scenario by SEM.

However, the proposed concept is general and can be adopted by other numerical methods such as FDM or FEM, and it can be directly applied to studying the hybrid numerical simulation of both the elastic wave equation and solid-fluid coupling, or Maxwell's equations.

Accepted Article

This article is protected by copyright. All rights reserved. where the first summation is assembled over all the surface elements e of the hybrid interface S, and the second one is over all the quadrature points x e p (p = 1,

is the total number of quadrature points on one surface element with the related quadrature weight α p . Subsequently, after introducing the discretization of Equation (A1) into the two surface integrals in Equation ( 8), the equivalent body force expression of equation (7) can be discretized as

In weak-form methods (such as the finite element method (FEM) or SEM), the traction source term in the weak form becomes

where w(x) is the abovementioned test function introduced in the second hybrid simulation. The moment-density tensor source term becomes

which is similar to the moment tensor implementation in [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF].

Appendix B Absorbing Boundary Condition implementation of the VM hybrid method

When the hybrid inputs are imposed on the local target model M l1 , the resultant potential q 1 on the hybrid interface S is

where q scatter is the scattered potential owing to the existence of heterogeneity in the local target model M l1 . The ABC [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF] can be used to absorb the scattered potential ∇q scatter on the hybrid interface S as follows:

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

where n is the unit outward normal of the boundary ∂Ω. With the hybrid inputs ∇q 0 and ∂ t q 0 , the second-order explicit Newmark time-marching scheme [START_REF] Newmark | A method of computation for structural dynamics[END_REF] can be used to perform the hybrid simulation in the local target model M l1 in the framework of SEM:

where the C 1 vector is obtained from the ABC, and Q1 and Q0 are the temporal derivatives of the potential vector calculated in the local target and reference models. M 1 is the global mass matrix, and K 1 is the global stiffness matrix assigned to the local target model M l1 . Note that Equation (B2) is used to absorb only the scattered waves and not the total wavefield. This boundary treatment naturally combines the forward wavefield with DSM or FK solutions and is the key part of the VM hybrid implementation.

Appendix C Best Scenario for the YM Method

In the best scenario, when the hybrid interface S completely follows the boundary of elements, namely, a Heaviside window function is adopted such that the diagonal value W of all the points in the local domain Ω l and on the hybrid interface S (e.g., the closed green line in Figure 2) is nonzero. The conceptual illustration of local target model M l1 associated with the YM method in this best scenario is depicted in Figure 1e, the blue part is where the hybrid input is imposed and the green part corresponds to the PML domain. Note that the blue part in Figure 1e for imposing the numerical hybrid inputs is outside the local target model M l1 . [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF] suggests that only part of values of the internal forces W e (K e 0 •Q e 0 ) and W e Q e 0 at the GLL nodes exactly on the hybrid interface S must be recorded according to the inner scheme in Equation ( 13) corresponding to the closed green line in the blue elements (hybrid elements in the BY hybrid method) in Figure 2. Because only the hybrid elements contribute to the numerical hybrid inputs, W e (K e 0 •Q e 0 ) and W e Q e 0 are related to only the GLL points on the hybrid interfaces under the selection of the element window function matrix W e . Note that in the best scenario, as in Figure 1e, although the blue domain is outside the hybrid interface S, it is an inner scheme as defined in [START_REF] Masson | On the numerical implementation of time-reversal mirrors for tomographic imaging[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved. The green waveform in (d1) is the enlarged sum of the blue and red dashed waveforms, and the green waveform in (d2) is the enlarged difference between the black waveform and the sum of the blue and red dashed waveforms.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth (a1, a2) is from the RP method. (b1, b2) is from the VM method. (c1, c2) is obtained from the MYM method. (d1, d2) is from the new method.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to JGR: Solid Earth

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. -53-

View publication stats