Acoustic wave-equation-based earthquake location

S U M M A R Y We present a novel earthquake location method using acoustic wave-equation-based traveltime inversion. The linear relationship between the location perturbation (δt0, δxs) and the resulting traveltime residual δt of a particular seismic phase, represented by the traveltime sensitivity kernel K(t0, xs) with respect to the earthquake location (t0, xs), is theoretically derived based on the adjoint method. Traveltime sensitivity kernel K(t0, xs) is formulated as a convolution between the forward and adjoint wavefields, which are calculated by numerically solving two acoustic wave equations. The advantage of this newly derived traveltime kernel is that it not only takes into account the earthquake–receiver geometry but also accurately honours the complexity of the velocity model. The earthquake location is obtained by solving a regularized least-squares problem. In 3-D realistic applications, it is computationally expensive to conduct full wave simulations. Therefore, we propose a 2.5-D approach which assumes the forward and adjoint wave simulations within a 2-D vertical plane passing through the earthquake and receiver. Various synthetic examples show the accuracy of this acoustic wave-equationbased earthquake location method. The accuracy and efficiency of the 2.5-D approach for 3-D earthquake location are further verified by its application to the 2004 Big Bear earthquake in Southern California.

. For Geiger's method, the relationship between the traveltime perturbation δt of a particular seismic phase (such as P-or S-wave) and the earthquake location perturbation (δt 0 , δx s ) can be written as where V xs is the velocity of the considered seismic phase at the source x s and l is the length parameter along the ray path (e.g. Thurber 1985;Ge 2003b;Chen et al. 2006;Thurber 2014). The Simplex method (e.g. Prugger & Gendzwill 1988;Rabinowitz 1988) is another primary algorithm for earthquake location (Ge 2003a,b). It utilizes a curve-fitting technique to search the minimums (optimal source locations) of the objective function in terms of the traveltime residuals. Different from Geiger's method, the Simplex method does not require derivative calculations (Prugger & Gendzwill 1988). It is also computationally efficient and stable (Rabinowitz 1988;Ge 2003b). More general reviews of various earthquake location methods, including the noniterative/iterative algorithms (Ge 2003a,b) or the methods for single-event/multiple-event location (Thurber 2014), can be found in Ge (2003a,b) and Thurber (2014). As discussed in many previous studies, the main factors determining the accuracy of earthquake location are the network geometry, measurement of phase arrival times, knowledge of background velocity model, and inversion algorithm (e.g. Waldhauser & Ellsworth 2000;Zhang et al. 2003;Schaff et al. 2004;Bondar et al. 2008;Maxwell et al. 2010). For the inversion algorithm, one of the key components is to accurately simulate seismic waves propagation (or at least predict the arrival times) in the given velocity model (e.g. Font et al. 2004;Chen et al. 2006). So far, most earthquake location methods are based on ray theory which assumes that traveltime and amplitude of an arrival only depend on (an)elastic properties along the geometrical ray path and ignores scattering, wave-front healing and other finite-frequency effects (e.g. Thurber 1985;Dahlen et al. 2000;Waldhauser & Ellsworth 2000;Zhang et al. 2003;Tong et al. 2011). However, seismic measurements such as traveltime and amplitude are sensitive to 3-D volumetric region off the ray path (e.g. Marquering et al. 1999;Dahlen et al. 2000;Tape et al. 2007). Furthermore, the ray theory is only valid when the scale length of the variation of material properties is much larger than the seismic wavelength (e.g. Rawlinson et al. 2010;Tong et al. 2014c,d). To accurately locate earthquakes, it is necessary to take into account the influence of off-ray structures by correctly capturing the interactions between the seismic waves and the potentially complex velocity model (e.g. Font et al. 2004;Liu et al. 2004). As pointed out by Tape et al. (2007) and Liu & Gu (2012), numerically solving the full wave equations is a promising approach to achieve this purpose. In the field of seismic tomography over the past decade, many high-resolution seismic images have been generated based upon numerical simulation of the full seismic wavefield, also known as full waveform inversion or adjoint tomogrpahy (e.g. Tape et al. 2009Tape et al. , 2010Fichtner & Trampert 2011;Zhu et al. 2012;Rickers et al. 2013;Tong et al. 2014d). Generally speaking, the main advantages of numerically solving the full wave equations in seismic inverse problems are the freedom to choose either a 1-D or 3-D velocity model and the accurate calculation of synthetic seismograms and sensitivity kernels for complex models which may help generate more reliable images (Kim et al. 2011;Liu & Gu 2012;Tong et al. 2014b).
Based upon 3-D numerical wave simulations with the spectral-element solvers (e.g. Komatitsch & Tromp 1999;Komatitsch et al. 2004), Liu et al. (2004) and Kim et al. (2011) developed and implemented a moment tensor inversion procedure to obtain focal mechanisms, depths, and moment magnitudes of earthquakes. In this study, we only focus on earthquake location and intend to incorporate the advantages of numerical solutions to the full wave equations. For the sake of simplicity, we use a simple acoustic wave equation to approximate the wavefield. The Fréchet derivative of traveltime with respect to earthquake location is derived based on a different procedure from Liu et al. (2004) and Kim et al. (2011). We call the proposed method the acoustic wave-equation-based earthquake location method.

THEORY
An earthquake that occurred at the source location x s is detected by a seismic station at x r . A perturbation δ m = (δt 0 , δx s ) in earthquake location m = (t 0 , x s ) generally leads to a traveltime shift δt of a particular seismic phase at the receiver x r . As a perturbation δt 0 in earthquake origin time t 0 directly causes an arrival-time shift δt 0 , only the traveltime shift δt − δt 0 caused by hypocentre perturbation δx s is considered in the following discussion. We assume that u(t, x r ) and s(t, x r ) are the synthetic seismograms calculated with earthquake location parameters m = (t 0 , x s ) and m + δ m = (t 0 , x s + δx s ), respectively. If |δm| |m|, we may expect that seismograms u(t, x r ) and s(t, x r ) are reasonably similar to each other. Based on this assumption, the traveltime shift δt − δt 0 of a particular phase can be calculated with high accuracy by maximizing the cross-correlation formula, where w(t) is a time window function for isolating the interested seismic phase in the time interval [0, T] (Tromp et al. 2005;Tong et al. 2014c). Based on the Born approximation (Dahlen et al. 2000), the maximum value of δt − δt 0 in eq. (2) can be alternatively computed via the following relationship Dashed and solid black curves denote the ray paths of the direct arrivals at x r generated at x s and x s + δx s , respectively. (b) An equivalent system in the velocity model c(x + δx s ) for the perturbated hypocentre x s + δx s and the receiver x r in panel (a). Identical seismic waves are generated at x s and recorded at x r − δx s . The black curve is the corresponding ray path of the direct arrival, same as the one in panel (a). where

Fréchet derivatives
We assume that seismic wave propagation with earthquake location parameter m = (t 0 , x s ) satisfies the acoustic wave equation where u(t, x) is the displacement field, c(x) represents the velocity distribution in the medium and f(t) denotes the source time function (if t ≤ 0 then f(t) = 0) at the source location x s . Similarly, displacement field s(t, x) = u(t, x) + δu(t, x) for an earthquake with location parameter m + δm = (t 0 , x s + δx s ) satisfies the following equation After a coordinate transformation from x to x + δx s (Alkhalifah 2010), it can be shown that seismic waves generated at hypocentre x s + δx s , propagated in the velocity model c(x) and recorded at receiver x r (Fig. 1a) are identical to the waves that are generated at x s , propagated in c(x + δx s ) and recorded at x r − δx s (Fig. 1b). The wavefield p(t, x) in velocity model c(x + δx s ) for a source at x s can be described as The relationship between s(t, x) and p(t, x) is Clearly, the only difference between eqs (6) and (4) lies in a velocity perturbation δx s · ∇c(x) from c(x) in eq. (4) to c(x + δx s ) in eq. (6), which accounts for the displacement perturbation δu(t, x) from u(t, x) to p(t, x) = u(t, x) + δu(t, x). Subtracting eq. (4) from eq. (6) and at University of Saskatchewan on March 1, 2016 http://gji.oxfordjournals.org/ Downloaded from neglecting the second-order terms, we obtain Multiplying an arbitrary test function q(t, x) on both sides of the first equation in eq. (8) and then integrating in the volume and the time interval [0, T], we have Integrating eq. (9) by parts gives If relationship (7) is substituted into eq. (3), then the traveltime shift δt − δt 0 can be expressed as Then by adding up eqs (10) and (11), using initial and boundary conditions in eq. (8), and assuming that the auxiliary field q(t, x) satisfies the following wave equation we obtain the relationship for traveltime shift as Define the Fréchet derivative G = G t 0 , G xs with respect to the origin time t 0 and hypocentre location x s as we can obtain the following relationship which links the traveltime shift δt and the earthquake location perturbation δm = (δt 0 , δx s ) = (δt 0 , δx, δy, δz) as The relationship (16) is similar to Geiger's method as shown in eq. (1), except that we need to solve two wave equations (4) and (12) to calculate the Fréchet derivatives in eq. (15) instead of tracing the ray path.

Inverse problem
For a particular seismic phase (such as direct P-wave or S-wave) generated by an event initially assumed at m = (t 0 , x s ), usually some time-shifts d = (δt i ) N × 1 can be measured at N seismic stations between the observed arrival times and the theoretical arrival times calculated in the known velocity model c(x). If the observing errors e = (e i ) N × 1 are taken into account, we have a concise matrix form for determining the earthquake location perturbation δm = (δt 0 , δx s ) as where G is an N × 4 sensitivity matrix with entries g ij . Here, g ij (j = 1, 2, 3, 4) are (G t 0 , G x , G y , G z ) in eq. (16) related to the ith observation. The standard damped least-squares solution for equation (17) with the minimum norm constraint is given by where I is a 4 × 4 identity matrix; λ is a non-negative damping parameter specified prior to the inversion to provide the intended weight of the minimum norm criterion. Once the earthquake location perturbation δ m is obtained, we can have an updated earthquake location m + δ m. Usually, an iterative procedure is needed to get an accuracy earthquake location. We call the whole process of calculating the Fréchet derivatives in eq. (16) and solving eq. (18) to obtain the final earthquake location 'the acoustic wave-equation-based earthquake location'.

Discussion on the Fréchet derivative G xs
The traveltime derivative G xs (x; x r , x s ) in eq. (15) with respect to the hypocentre position x s consists of two terms G 1 xs and G 2 xs as and Generally speaking, G 1 xs mainly accounts for the influence caused by the geometrical mislocation of the hypocentre (zeroth-order influence), while G 2 xs considers the impact of the spatial variation of the velocity model (first-order influence). For example, if the hypocentre moves closer to the receiver, G 1 xs indicates that the traveltime should be reduced. But even if the hypocentre is closer to the receiver, it is possible that the considered seismic phase may pass through a region with a lower velocity structure. In this case, G 2 xs suggests an increase in traveltime. Specifically, we consider a homogeneous model with constant velocity c. Assuming that an impulsive source f(t) = δ(t − t 0 ) is exerted at the source, the wavefield recorded at x r has an analytical form as u(t, , and the first term in eq. (15) is simplified to G 1 xs = −(x r − x s )/(c|x r − x s |) and the second term G 2 xs vanishes. G 1 xs is the same as the derivative of arrival time with respect to the hypocentre location as shown in eq. (1), which is derived in the framework of ray theory (Geiger 1910;Engdahl & Lee 1976;Thurber 2014). That is to say, G 1 xs can be roughly considered as the Fréchet derivatives described in Geiger's method for a homogeneous model; and the other term G 2 xs additionally takes into account the velocity variation, which could be helpful for earthquake location in complex media. Sensitivity or Fréchet kernels defined as the volumetric densities of the Fréchet derivative are widely used in tomographic studies of volumetric material properties (e.g. Dahlen et al. 2000;Tromp et al. 2005;Fichtner & Trampert 2011;Tong et al. 2014a). They have direct physical meanings in interpreting the finite-frequency effects of propagating seismic waves. Based on G 2 xs in eq. (20), we can also introduce sensitivity or Fréchet kernel for earthquake location as The velocity model c(x) and its spatial variation ∇c(x), the distribution of the hypocentre x s and receiver x r , and the considered seismic phase together determine the sensitivity kernel (21) in earthquake location inversion. Note that the Fréchet kernel for earthquake location is simply multiplying ∇c(x) to the volumetric sensitivity kernel for c(x) in tomographic inversions (Tong et al. 2014d).

2-D forward modelling and 3-D inversion
It is known that seismic tomography based on the simulation of full wavefield is computationally expensive (e.g. Chen et al. 2007;Tape et al. 2009;Zhu et al. 2012;Tong et al. 2014d). Since two wave equations need be solved to obtain the Fréchet derivative, the acoustic wave-equation-based earthquake location method developed in this study is also computationally demanding, especially when 3-D spatial geometry is considered. To reduce the computational cost, following the approximation technique used in wave-equation-based traveltime seismic tomography (Tong et al. 2014c,d), we suggest to restrict the simulation of wavefield in a 2-D vertical plane passing through the hypocentre and receiver. The horizontal direction from the hypocentre to receiver is defined as r. The angle between r and the x-axis is θ . Within the 2-D vertical plane, the traveltime Fréchet derivative with respect to the hypocentre location x s has two components G xs = (G r , G z ). To invert for 3-D earthquake location, we can further project the component G r onto two horizontal directions as (G x , G y ) = (G r cos θ, G r sin θ ). In this scenario, the linearized relationship between traveltime perturbation and earthquake location perturbation shown in eq. (16) can be modified into Eq. (22) only requires the simulation of full wavefield in a 2-D vertical plane but can be used to invert for 3-D earthquake location. This 2.5-D strategy is computationally more efficient compared to the one based on 3-D forward modelling.

N U M E R I C A L E X A M P L E S
We test the performances of the acoustic wave-equation-based earthquake location method with four synthetic numerical examples.

The two-layer model
We first consider a two-layer velocity model (Fig. 2a). For demonstration purpose, only 2-D geometries are considered in this and the next proof-of-concept examples. The implementation of the proposed earthquake location method in 3-D geometries is very similar. As an attempt to fully demonstrate the performances of the new earthquake location method, three different cases are discussed: (1) The initial hypocentre location is in the top layer, but the true location is in the bottom layer.
(2) Both initial and actual hypocentre locations are in the bottom layer.
(3) The initial location is in the bottom layer, while the true location is in the top layer (Fig. 2a). The origin time is not perturbed but iteratively updated in this example. Ten seismic stations with an equal spacing of 1.0 km are employed at the surface to record seismograms generated Table 3. The same as Table 1 but for the third case in Fig. 2(a).   at the initial, iteratively updated and actual earthquake locations. In this study, only the first P-wave arrivals are used to locate the earthquakes and all the synthetic seismograms are generated with a Ricker wavelet point source f(t − t 0 ) in eq. (4) which has an analytical form of f 0 is the dominant frequency chosen as f 0 = 1.0 Hz here and A is the normalization factor. A high-order finite-difference method is adopted to simulate the full wavefields (Tong et al. 2014c). Fig. 2(a) and Tables 1-3 show the iteratively updated earthquake locations. The corresponding traveltime residuals of the first P-arrivals at each station are shown in Figs 2(b)-(d). We can observe that the acoustic wave-equation-based earthquake location method can accurately locate the three earthquakes after about 4 iterations, showing the validity of this new technique. To examine the effect of velocity variation on earthquake location, we calculate and show one Fréchet kernel K xs in Fig. 3. Since the velocity model has no lateral variation, the horizontal component of G 2 xs is zero. Meanwhile, the vertical component of G 2 xs reflects the influence of the velocity discontinuity on locating the depth of the earthquake.

Inversion with seismic data of different frequencies
As discussed previously, the acoustic wave-equation-based earthquake location method naturally takes into account the finite-frequency effects of the seismic waves. In the second example, we explore the influence of different frequency contents of seismic waves on the accuracy of earthquake location in a 3-layer velocity model with lateral variation (Fig. 4). The hypocentre is actually located in the top layer but initially assumed to be near the bottom of the middle layer (Fig. 4). We use the first P-waves at dominant frequencies f 0 = 5.0 Hz and f 0 = 1.0 Hz to locate this earthquake separately. Fig. 4 shows the iteratively updated hypocentre locations and the corresponding root mean square (RMS) of the first P-wave traveltime residuals. Tables 4 and 5 quantitatively demonstrate the earthquake location throughout the iterations for using the high-frequency and low-frequency data, respectively. Generally speaking, using both high-and low-frequency seismic data can accurately locate the earthquake with the acoustic wave-equation-based earthquake location method. But a closer examination reveals  that higher-frequency seismic data can locate the earthquake more precisely. For example, after 12 iterations the RMS of the first P-wave traveltime residuals is in the order of 10 −4 s (Fig. 4c) and 10 −3 s (Fig. 4d) for the high-and low-frequency cases, respectively. Tables 4 and 5 also indicate that higher-frequency data can find a closer earthquake location to the 'true' location after the same number of iterations. Fig. 5 displays the sensitivity kernels for the horizontal location and depth, which are corresponding to a single source-receiver pair. The obvious observation is that the high-frequency data have a smaller influence zone (sensitivity kernel with nontrivial value) than the low-frequency data do. That is to say, different frequency components of the same seismic phase may have slightly different arriving times if the velocity model has spatial variation. This also suggests that we need to consider both the finite-frequency effects of the propagating seismic waves and the velocity variation for accurate earthquake locations as the proposed acoustic wave-equation-based technique does.

Earthquake location with the 2.5-D approach
Previous numerical examples have verified the concept of acoustic wave-equation-based earthquake location in 2-D geometries. In Section 2, we also proposed to use a 2.5-D approach to locate earthquakes in 3-D space with 2-D forward modellings to reduce the computational cost. The third example is designed to test this 2.5-D approach in an extreme case. The considered velocity model consists of the crust and the mantle, containing an undulated Moho and a subduction zone with a thin low velocity layer atop a fast velocity layer (Fig. 6). The central part of the Moho has a maximum elevation of 10.0 km from the flat position at the depth of 30.0 km. Comparing to the surrounding mantle, the velocities of the fast and slow velocity layers of the subduction zone are perturbed by +4 per cent and −6 per cent, respectively. The initial hypocentre location of the earthquake is in the high velocity layer of the subduction zone, about 115 km away from the 'true' location in the above low velocity zone (Fig. 6). The origin time of the earthquake is also delayed by 5.0 s from the 'actual' origin time. Three arrays of and a total of 21 seismic stations on the surface are used to record seismograms (Fig. 6). The dominant frequency of all the seismograms is  Table 6. The iteratively updated earthquake location obtained by using the 2.5-D approach in a crust-over-mantle model with an undulated Moho and a subduction zone (Fig. 6). f 0 = 0.2 Hz. Fig. 6 displays the iteratively relocated hypocentre locations (Figs 6a-c) and the corresponding root mean square values of the first P-wave traveltime residuals (Fig. 6d). Meanwhile, Table 6 quantitatively demonstrates the initial, updated, and target earthquake locations. The hypocentre location and origin time errors at the 20th iteration are less than 0.06 km and 0.005 s, respectively. The errors are mainly caused by the comparatively large grid interval (1.0 km) used in the forward modelling scheme, the numerical wave-equation solver itself, and the resolving ability of the relatively low frequency (f 0 = 0.2 Hz) data. All the results of this example indicate that the 2.5-D version of the acoustic wave-equation-based earthquake location method can accurately locate the earthquake even though the initial location is far away from the 'true' location.

Influence of inaccurate velocity model
The accuracy of earthquake location relies on how well we know the velocity model. Since the 'exact' velocity model is always not available, we should bear in mind that earthquake locations are usually obtained by using approximate velocity models in real applications. All the previous examples reveal that the acoustic wave-equation-based earthquake location method has excellent performance in locating earthquakes if the known velocity models are exactly the ones in which the data are generated. To investigate our technique for locating real earthquakes, in the fourth example we test the performance of the 2.5-D version of the proposed method when only an approximate model is provided. For a direct comparison purpose, we relocate the earthquake discussed in the third example again. All the parameters are the same except that the velocity structure for generating synthetic seismograms is simplified into a two-layer model with a flat Moho at the depth of 30.0 km (Fig. 7). Note that the data generated at the 'true' earthquake location is still calculated in the 'true' velocity model as that of Figs 6(a)-(c). The iteratively updated earthquake locations relocated in the simple two-layer velocity model are shown in Fig. 7 and Table 7. After 20 iterations, the inverted earthquake location is close to the 'actual' location but there is still an obvious gap between them (Fig. 7). Table 7 reveals that the final spatial location error is about 10.0 km, though it is much smaller than the initial location error of 115.0 km. Specifically, the earthquake depth is 7.73 km deeper than the real depth and the origin time t 0 is advanced by 0.48 s in contrast to initially being delayed by 5.0 s. Comparing the 'true' and 'approximate' velocity models suggests that the velocity model used for earthquake location does not account for the crustal low velocity structure below 30.0 km and the low velocity layer of the subducting slab. The failure of considering these low velocity structures makes the earthquake deeper. This example indicates that relocating earthquakes in approximate velocity models is likely to obtain approximate earthquake locations with the acoustic wave-equation-based earthquake location method.

R E A L D ATA A P P L I C AT I O N
The main purpose of developing the acoustic wave-equation-based earthquake location method is to use it to locate earthquakes in real applications. Liu et al. (2004) proposed a moment tensor inversion procedure based upon spectral-element simulations and obtained focal mechanisms, depths, and moment magnitudes of three southern California earthquakes including the 2004 M L 5.4 Big Bear event (Fig 8a). They claimed that their inversion results are generally in good agreement with estimates based upon traditional body-wave and surface-wave inversions . In this study, we also choose the 2004 Big Bear earthquake as the test event to further explore our new 2.5-D earthquake location method. Table 7. The same as Table 6 but the results are obtained in a simple crust-over-mantle model (Fig. 7).  We download broadband seismic data from Southern California Seismic Network (http://scsn.org). Only the first P-waves recorded at the 10 nearest seismic stations are used to locate the 2004 Big Bear earthquake. All the 10 seismic stations are at epicentral distances less than 45 km (Fig. 8a). The earthquake location is conducted in a 1-D layered model, which is separated by discontinuities at the depths of 5.5, 16.0 and 29.2 km (Tong et al. 2014d). The P-wave velocities in the four layers are 5.5, 6.3, 6.7 and 7.8 km s −1 , respectively. The dominant frequency of the source time function (23) is chosen as f 0 = 1.0 Hz. As a result, the nontrivial part of each synthetic seismogram lasts for nearly 2.0 s (twice the dominant period). To be consistent, 2-s time windows starting at the picked onset times of the P-arrivals (from the SCSN catalogue) are used to isolate the portions of the data for earthquake location. A Butterworth filter between 0.1 and 1.3 Hz is then applied to the windowed seismograms. For effectively measuring the cross-correlation traveltime residuals between the synthetics and observed data, all the 2-D waveforms are converted to 3-D seismograms using the conversion formula derived in Miksat et al. (2008). We relocate the Big Bear earthquake starting from its SCSN catalogue location. Similar to Liu et al. (2004), the origin time of the earthquake is fixed and only the hypocentre location is updated in this test. Table 8 contains the iteratively inverted hypocentre locations and the corresponding RMS values of the first P-wave traveltime residuals. We can observe that the lateral location (Longitude, Latitude) of the Big Bear earthquake has very small variations but the depth becomes deeper by more than 5.0 km throughout a total of 14 iterations. The RMS value of the traveltime residuals is reduced from 0.33 to 0.15 s with a minimum value of 0.13 s at the first iteration. Although the RMS value reaches a minimum after a single iteration, we do not stop the relocation procedure because the earthquake location still varies significantly in the following iterations. After 14 iterations, we almost obtain a stable solution and consider it as the final earthquake location of the Big Bear event. Using three different inversion schemes including the spectral-element moment tensor inversion, Liu et al. (2004) concluded that the depth of the 2004 Big Bear earthquake is about 6.4 ± 0.2 km. The depth obtained by the acoustic wave-equation-based earthquake location method is 6.35 km (Table 8), consistent with the result based upon the spectral-element moment tensor inversion. Fig. 9 shows the results of waveform fitting at two representative stations. We can observe that the synthetic waveforms generated at the final earthquake location have smaller phase shifts from the recorded data than the ones corresponding to the initial location do. This implies that we have obtained a more accurate location for the Big Bear earthquake. In addition, we can investigate the relationship between earthquake occurrence and seismic velocity providing that the earthquake is accurately located. Tong et al. (2014d) mapped the seismic velocity structures of the 1992 Landers earthquake area and the 2004 Big Bear earthquake is actually in the investigated region (Fig 8a). A vertical view of the P-wave velocity model recovered by Tong et al. (2014d) which passes through the 2004 Big Bear earthquake and the 1992 Landers earthquake is shown in Fig. 8(b). We can observe that the Big Bear earthquake is initially located in a high Vp anomaly close to the surface but its final location is in a transition zone between high Vp and low Vp structures. The inverted earthquake location with the acoustic wave-equation-based method is in agreement with the statement that many large or relatively large crustal earthquakes occurred in regions with significant seismic property variations (Tong et al. 2014d). The tomographic evidence also indicates the validity of the newly proposed earthquake location method.

D I S C U S S I O N A N D C O N C L U S I O N S
The full wavefield simulation honours the finite-frequency effects of the propagating seismic waves and accurately captures the interactions between seismic waves and complex (such as strongly heterogeneous) structures. These advantages prompted us to derive the acoustic waveequation-based earthquake location method. Considering that the full 3-D seismic numerical modelling is always computationally demanding or even prohibitive, a 2.5-D version of this method is designed to lessen the computation burden. From the synthetic examples in Section 3, we (1) The newly proposed method can precisely locate earthquake location if high-quality data and exact velocity model are ideally available.
(2) High-frequency data has better resolving ability than low-frequency data. (3) If only an approximate velocity model is known, the acoustic wave-equation-based method probably finds a location which is not but around its actual location. Furthermore, the validity of the acoustic wave-equation-based earthquake location method in real applications is verified by applying it to relocate the 2004 Big Bear earthquake and comparing the result with the ones generated by other techniques such as the spectral-element moment tensor inversion of Liu et al. (2004). The results of this study suggest that the 2.5-D version of the acoustic wave-equation-based earthquake location method can be independently used to locate earthquakes with computational efficiency and accuracy. However, it still needs more computational resources than Geiger's method does. As discussed in Section 2, the main advantage of the acoustic wave-equation-based technique is that it can accurately calculate the Fréchet derivatives, especially in complex velocity models. Since the earthquake location is iteratively updated, in future studies we can first use Geiger's method to get an initial location and then adopt the acoustic wave-equation-based method to refine the location. The combination of these two methods should reduce the computation cost. Furthermore, the velocity model is a key component in earthquake location. The extent to which we know the velocity model determines how accurately we can locate the earthquake. For a sole earthquake location inversion, the error in velocity model may trade-off with the inverted results. Therefore, a simultaneous inversion for velocity model and earthquake location may be necessary for obtaining a more accurate earthquake location. In addition, the 2.5-D approach suggests using 2-D forward modelling for 3-D hypocentre location. If the velocity model is very complex and the off-plane effects cannot be ignored, we can resort to the acoustic wave-equation-based earthquake location method based upon 3-D forward modelling. For the sake of conciseness, we put no emphasis on location uncertainty in this study. Readers can refer to Thurber (2014) for a thorough discussion on how to analyse the precision and accuracy of earthquake location.
In conclusion, the proposed acoustic wave-equation-based earthquake location method provides a new efficient way for accurately locating earthquakes especially when complex velocity models are present. We expect that it can be used in seismic tomography studies with passive data, earthquake source parameters (such as centroid-moment tensor) inversion (Kim et al. 2011), microseismic monitoring in the process of hydraulic fracturing, nuclear explosion monitoring, and many other fields. For the next stage, this acoustic wave-equationbased earthquake location technique can be further developed to take advantage of the full waveform content instead of using only the cross-correlation traveltime information as in this study.  Wessel & Smith 1991). We thank Egill Hauksson and two anonymous reviewers for providing constructive comments and suggestions that improved the manuscript.