Seismic Traveltime Tomography Research Articles

Application of travel-time seismic tomography for archaeological studies of building foundations and basements

The applicability of the travel-time seismic tomography algorithm with the adaptive environment parameterization for the engineering and archaeological studies is investigated. The investigation object is the environment under the ruins of the St. Onuphrius Church at the Big Solovetsky Island. Experimental geometry for the tomography studies is designed and the data are acquired and processed. The goal is to study in detail the church basement and underlying geological environment structure that can not be recovered with other geophysical methods. The synthetic experiments are performed using the standard checkerboard tests with the cells of 2, 3, and 4 m in size in order to estimate the adequacy of the results obtained. The recovered three-dimensional velocity pattern is in a good agreement with the data on the church and archaeological information.

Crustal P-wave velocity structure of the Longmenshan region and its tectonic implications for the 2008 Wenchuan earthquake

The P-wave velocity structure of the crust in the Longmenshan region has been imaged by seismic travel time tomography using local and regional first P-wave arrivals recorded from 2000 to 2008. The tomographic model provides a way to analyze the deep tectonics of the Longmenshan fault belt and the tectonic implications for the 2008 Ms8.0 Wenchuan earthquake. The P-wave velocity images indicate that the initial rupture site and focal depth of the Wenchuan earthquake, together with the direction of rupture propagation, closely relate to the crustal structure of the Longmenshan region. The Pengguan massif to the west of the Longmenshan fault belt is characterized by high velocity anomalies, suggesting that the crust has a strong strain strength that can accumulate large stresses over a long period. The Ms8.0 Wenchuan earthquake is located at the southwestern end of the Pengguan massif and the western edge of the Sichuan Basin. The collision between the Pengguan massif and the Sichuan Basin becomes the primary reason for the occurrence of the Ms8.0 Wenchuan earthquake. To the north of Wenchuan, the occurrence and propagation of rupture benefit from low velocity anomalies along the Longmenshan fault belt; whereas to the south of Wenchuan, the brittle rupture can occur with more difficulty in relatively weak crust with low velocities. This may be one of the reasons for the absence of aftershocks to the south of Wenchuan, and the rupture induced by the Ms8.0 Wenchuan earthquake propagating from the north to the south along the Longmenshan fault belt. The deep geodynamics of the Ms8.0 Wenchuan earthquake may occur due to the discrepancy of crustal structures on the two sides of the Longmenshan fault belt. Ductile deformation and crustal flow can easily occur in the weak middle-lower crust beneath the Songpan-Garze orogenic belt. The eastward movement of the Tibetan Plateau is obstructed by the rigid lithosphere of the Sichuan Basin, and then the thickening of the middle-lower crust and vertical deformation occur in the crust of the Longmenshan fault belt. In addition, the down-warping of the Moho and the basement thrusting onto the range front induced crustal deformation and strain accumulation, which provided the potential energy to trigger the occurrence of the Ms8.0 Wenchuan earthquake.

Traveltime seismic tomography with adaptive wavelet parameterization

An algorithm for solving the inverse kinematic problem of traveltime seismic tomography is developed and tested. The algorithm is intended for imaging the three-dimensional (3D) velocity model composed of a layer underlain by a half-space. This algorithm considers the bottom boundary of the layer as a first-order seismic velocity discontinuity with unknown position that has to be determined in the inversion together with the velocity variations inside the overlying layer and the sub-interface boundary velocities. The inversion can be applied to the travel times of refracted, head and reflected waves. The main idea behind the algorithm is the adaptive parameterization of the medium by the sparse Haar wavelet series expansion. In order to throw off the poorly resolved coefficients of expansion, we suggest using two empirical local resolution measures: the number of seismic rays crossing the support of the corresponding wavelet support area and their angular coverage, i.e., the spread in the azimuths of these rays. The adequacy of these measures is tested by their comparison with the estimation of the diagonal elements of the resolution matrix on the synthetic examples. This comparison proved that the proposed measures can be successfully applied for statistical estimation of the resolution and for constructing the adaptive parameterization. It was shown also that the best results are achieved while using the number of rays normalized to the size of the wavelet support together with their angular coverage. An automated procedure for throwing off poorly resolved unknowns is developed. The parameters of this procedure can be tuned to provide the desired level of detail of the model to be reconstructed. The synthetic checkerboard testing proved the efficiency of the algorithm. The proposed algorithm can be applied to solve different types of problems, including regional seismic studies, as well as exploration and engineering seismology. The use of this algorithm is especially convenient when the medium is essentially three-dimensional and when the conventional seismic methods implying regular network measurements directly above the studied structure (such as the common depth point method) are inapplicable, e.g., in the seismic studies of the foundations of buildings and in rugged terrains.

Recent advances in optimized geophysical survey design

Survey design ultimately dictates the quality of subsurface information provided by practical implementations of geophysical methods. It is therefore critical to design experimental procedures that cost effectively produce those data that maximize the desired information. This review cites recent advances in statistical experimental design techniques applied in the earth sciences. Examples from geoelectrical, crosswell and surface seismic, and microseismic monitoring methods are included. Using overdetermined 1D and 2D geoelectrical examples, a minor subset of judiciously chosen measurements provides a large percentage of the information content theoretically offered by the geoelectrical method. In contrast, an underdetermined 2D seismic traveltime tomography design study indicates that the information content increases almost linearly with the amount of traveltime data (source-receiver pairs) considered until the underdeterminancy is reduced substantially. An experimental design study of frequency-domain seismic-waveform inversion experiments reveals that a few optimally chosen frequencies offer as much subsurface information as the full bandwidth. A nonlinear experimental design for a seismic amplitude-versus-angle survey identifies those incidence angles most important for characterizing a reservoir. A nonlinear design example shows that designing microseismic monitoring surveys based on array aperture is a poor strategy that almost certainly leads to suboptimal designs.

On Mathematical Aspects of a Combined Inversion of Gravity and Normal Mode Variations by a Spline Method

This paper provides a brief overview of two linear inverse problems concerned with the determination of the Earth’s interior: inverse gravimetry and normal mode tomography. Moreover, a vector spline method is proposed for a combined solution of both problems. This method uses localised basis functions, which are based on reproducing kernels, and is related to approaches which have been successfully applied to the inverse gravimetric problem and the seismic traveltime tomography separately.

Regularization of linear and non-linear geophysical ill-posed problems with joint sparsity constraints

Summary In this paper, we deal with the solution of linear and non-linear geophysical ill-posed problems by requiring the solution to have sparse representations in two appropriate transformation domains, simultaneously. Geological structures are often smooth in properties away from sharp discontinuities (i.e. jumps in 1-D and edges in 2-D). Thus, an appropriate ‘regularizer’ function should be constructed so that recovers the smooth parts as well as the sharp discontinuities. Sparsity inversion techniques which require the solution to have a sparse representation with respect to a pre-selected basis or frames (e.g. wavelets), can recover the regions of smooth behaviour in model parameters well, but the solution suffers from the pseudo-Gibbs phenomenon, and is smoothed around discontinuities. On the other hand, requiring sparsity in Haar or finite-difference (FD) domain can lead to a solution without generating smoothed edges and the pseudo-Gibbs phenomenon. Here, we set up a regularizer function which can be benefited from the advantages of both wavelets and Haar/FD operators in representation of the solution. The idea allows capturing local structures with different smoothness in the model parameters and recovering smooth/constant pieces of the solution together with discontinuities. We also set up an information function without requiring the true model for selecting optimum wavelet and parameter β which controls the weight of the two sparsifying operators in the inverse algorithm. For both linear and non-linear geophysical inverse problems, the performance of the method is illustrated with 1-D and 2-D synthetic examples and a field example from seismic traveltime tomography. In all of the examples tested, the proposed algorithm successfully estimated more credible and high-resolution models of the subsurface compared to those of the smooth and traditional sparse reconstructions.

Modified Finslerian spaces for more effective ray-tracing. Application to seismic travel-time tomography

Abstract This paper deals with a new ray-tracing approach for seismic travel-time tomography, using the Finslerian geometry. Traditional ray-tracing methods produce poor results because non-uniform coverage of rays in the whole computational domain. We developed and tested the new algorithm based on the Finsler geometry. Tests were carried out using both synthetic and real data. The results are described in this paper.

3D density model of the Central Andes

We developed a 3D density model of the continental crust, the subducted plate and the upper mantle of the Central Andes between 20–29°S and 74–61°W through the forward modelling of Bouguer anomaly. The goal of this contribution is to gain insight on the lithospheric structure integrating the available information (geophysical, geologic, petrologic, and geochemical) in a single model. The geometry of our model is defined and constrained by hypocentre location, reflection and refraction on and offshore seismic lines, travel time and attenuation tomography, receiver function analysis, magnetotelluric studies, thermal models and balanced structural cross-sections. The densities allocated to the different bodies are calculated considering petrologic and geochemical data and pressure and temperature conditions. The model consists of 31 parallel E–W vertical planes, where the continental crust comprises distinct bodies, which represent the different morphotectonic units of the Central Andes. We include a partial melting zone at midcrustal depths under the Altiplano-Puna (low-velocity zone) and consider the presence of a rheologically strong block beneath the Salar de Atacama basin, according to recent seismic studies. Contour maps of the depth of the continental Moho, the thickness of the lower crust and the depth to the bottom of the lithosphere below South America are produced. The possible percentage of partial melt in the Central Andes low-velocity zone is estimated. The residual anomaly is calculated by subtracting from the Bouguer anomaly the gravimetric effect of the modelled subducted slab and of the modelled Moho. Isostatic anomalies are calculated from regional and local isostatic Mohos calculated with and without internal loads, derived from our gravity model, which are then compared to the modelled continental Moho. This study contributes to a more detailed knowledge of the lithospheric structure of this region of the Andes and provides an integrated 3D density model, which may be used by the geoscientific community as a tool that can help to understand and interpret the geodynamic features and processes acting along the Central Andes.

Near-surface seismic traveltime tomography using a direct-push source and surface-planted geophones

Information about seismic velocity distribution in heterogeneous near-surface sedimentary deposits is essential for a variety of environmental and engineering geophysical applications. We have evaluated the suitability of the minimally invasive direct-push technology for near-surface seismic traveltime tomography. Geophones placed at the surface and a seismic source installed temporarily in the subsurface by direct-push technology quickly acquire reversed multioffset vertical seismic profiles (VSPs). The first-arrival traveltimes of these data were used to reconstruct the 2D seismic velocity distribution tomographically. After testing this approach on synthetic data, we applied it to field data collected over alluvial deposits in a former river floodplain. The resulting velocity model contains information about high- and low-velocity anomalies and offers a significantly deeper penetration depth than conventional refraction tomography using surface-planted sources and receivers at the investigated site. A combination of refraction seismic and direct-push data increases resolution capabilities in the unsaturated zone and enables reliable reconstruction of velocity variations in near-surface unconsolidated sediments. The final velocity model structurally matches the results of cone-penetration tests and natural gamma-radiation data acquired along the profile. The suitability of multiple rapidly acquired reverse VSP surveys for 2D tomographic velocity imaging of near-surface unconsolidated sediments was explored.

Bayesian analysis of data and model uncertainty in 3D seismic travel-time tomography

Abstract We present a probabilistic analysis of seismic travel-time equations using the Bayesian Method. The assessment of models and data is crucial in 3D seismic travel-time tomography, and a method quantitatively assess the quality of both the data and the model is necessary in order to attain the most realistic results. The Bayesian method that we propose here is more effective than the frequentist approach, both in analysis time and uncertainty minimization, when processing large sets of tomographic data.

A continuation approach to regularization of ill-posed problems with application to crosswell-traveltime tomography

In most geometries in which seismic-traveltime tomography is applied (e.g., crosswell, surface-reflection, and VSP), determination of the slowness field using only traveltimes is not a well-conditioned problem. Nonuniqueness is common. Even when the slowness field is uniquely determined, small changes in measured traveltimes can cause large errors in the computed slowness field. A priori information often is available — well logs, initial rough estimates of slowness from structural geology, etc. — and can be incorporated into a traveltime-inversion algorithm by using penalty terms. To further regularize the problem, smoothing constraints also can be incorporated using penalty terms by penalizing derivatives of the slowness field. What weights to use on the penalty terms is a major decision, particularly the smoothing-penalty weights. We use a continuation approach in selecting the smoothing-penalty weights. Instead of using fixed smoothing-penalty weights, we decrease them step by step, using the slowness model computed with the previous, larger weights as the initial slowness model for the next step with the new, smaller weights. This continuation approach can solve synthetic problems more accurately than does one that uses fixed smoothing-penalty weights, and it appears to yield more features of interest in real-data applications of traveltime tomography. We have formulated guidelines for making the many choices needed to implement this continuation strategy effectively and have developed specific choices for crosswell-traveltime tomography.

P wave seismic velocity and Vp/Vs ratio beneath the Italian peninsula from local earthquake tomography

We investigate the P wave velocity structure and the Vp/Vs ratio beneath the Italian peninsula down to 100 and 60 km depth respectively by seismic travel time tomography. We invert data provided by the International Seismological Centre (ISC) (1997–2005), making use of some alternative strategies for the travel time approach in a well constrained and worldwide adopted code (SIMULPS). Resolution for the different layers is discussed and sensitivity analyses are performed through test inversions to explore the resolution characteristics of the model at different spatial scales. The resulting tomographic images provide a detailed sketch of the P wave anomalies, clearly showing, among the other features, the shape of the Ivrea body in the Western Alps, the upwelling of the oceanic crust in the Ligurian Sea and the slab under the Calabrian arc. They are less informative for the Vp/Vs ratio. Nevertheless, some features are very interesting and deserve further investigation like the anomalous decrease of the Vp/Vs ratio under the Ligurian Sea or the variations of the Vp/Vs ratio calculated in the first 10 km depth of the Apenninic region with respect to the lower values of the Alpine region at the same depth. The tomographic cross sections reveal a continuous superposition of two kinds of crusts (transitional over Adriatic) all along the peninsula but do not show any slab, intended as a clear, vertical downgoing high velocity material in either the northern or central Apennines.

High-resolution seismic characterization in an urban area: Subway tunnel construction in Barcelona, Spain

Traveltime seismic tomography used to aid subway tunnel drilling in Barcelona, Spain, provides a detailed characterization of the shallow subsurface, including a complex network of faults and dikes. We encountered many difficulties while working in this urban setting, such as cultural noise and the inability to deploy instruments along much of the tunnel trace because of the street layout. Furthermore, the shallow subsurface is highly heterogeneous (because of gas and water lines, building foundations, and the sewage system), obscuring deeper structures. We acquired seismic data using a 10-s vibroseis sweep that provided relatively good signal-to-noise (S/N) ratio. First-arrival traveltime tomography proved to be a key tool to differentiate rock vol-umes characterized by different physical properties — in particular, the contact between the weathered upper layer and the more competent rock below. We interpreted the seismic results using surface geologic observations and borehole cores. Low seismic velocities [Formula: see text] characterize the shallow part of the subsurface, representing quaternary sediments and weathered granite. Several high-velocity anomalies (up to [Formula: see text]) are interpreted as porphyritic dikes and are surrounded by low-velocity anomalies that correspond with fault systems that cut and displace the dikes. The tomographic velocity models also provide a detailed image of an important Miocene fault that intersects the tunnel trace. These results provide insights important for the drilling process.

3D seismic traveltime tomography imaging of the shallow subsurface at the C O2 SINK project site, Ketzin, Germany

A 3D reflection seismic survey was performed in 2005 at the Ketzin carbon dioxide [Formula: see text] pilot geological-storage site (the [Formula: see text] project) near Berlin, Germany, to image the geological structure of the site to depths of about [Formula: see text]. Because of the acquisition geometry, frequency limitations of the source, and artefacts of the data processing, detailed structures shallower than about [Formula: see text] were unclear. To obtain structural images of the shallow subsurface, we applied 3D traveltime tomography to data near the top of the Ketzin anticline, where faulting is present. Understanding the shallow subsurface structure is important for long-term monitoring aspects of the project after [Formula: see text] has been injected into a saline aquifer at about [Formula: see text] depth. We used a 3D traveltime tomography algorithm based on a combination ofsolving for 3D velocity structure and static corrections in the inversion process to account for artefacts in the velocity structure because of smearing effects from the unconsolidated cover. The resulting velocity model shows low velocities of [Formula: see text] in the uppermost shallow subsurface of the study area. The velocity reaches about [Formula: see text] at a depth of [Formula: see text]. This coincides approximately with the boundary between Quaternary units, which contain the near-surface freshwater reservoir and the Tertiary clay aquitard. Correlation of tomographic images with a similarity attribute slice at [Formula: see text] (about [Formula: see text] depth) indicates that at least one east-west striking fault zone observed in the reflection data might extend into the Tertiary unit. The more detailed images of the shallow subsurface from this study provided valuable information on this potentially risky area.

Using seismic traveltime tomography in geoarchaeological exploration: an application at the site of Chatby cemeteries in Alexandria, Egypt

We apply seismic traveltime tomography in an investigation of the subsurface of a geoarchaeological site at Chatby cemeteries (Alexandria, Egypt). We are interested in obtaining information about the geological structure and the remains of man‐made structures. Energy sources are placed at the surface and recorders are distributed vertically in an underground location. Traveltime data are inverted with a three‐dimensional tomographic algorithm. Inversion results are evaluated and interpreted with the aid of independent geophysical and archaeological information, and various tests are carried out with different initial velocity models and synthetic traveltimes.The inversion results show good correlation with the estimated geological structure, in particular at deeper levels, as suggested by independent geophysical and archaeological information and later confirmed by excavation. The depth and variations of the bedrock surface are satisfactorily imaged. Indications of possible man‐made structures are presented and discussed. Shallow velocity features are evaluated and interpreted using tests with synthetic traveltimes. The main goal of imaging the deeper structure was successful.As inversion tests with different initial velocity models indicate, the subsurface structure and the complexity of data information in this particular case are better handled when using a single‐layer model than when using an average two‐layer model. The single‐layer model allows for a more versatile velocity variation and raypath distribution, while the two‐layer model leads to an inversion result that is considerably influenced by the initial velocity distribution. Enhancement of particular aspects of the single‐layer inversion result can be derived with the use of an initial velocity model, that permits a gradual increase of velocity with depth. The gradient model improves the raypath distribution and limits velocity artefacts.The resulting images can assist in the assessment of the archaeological potential of a site, particularly in excavation design, in a non‐destructive manner. They can be also of great value in other applications, for example in engineering works, when detailed information on the subsurface is sought in areas of restricted access.

Small-scale convection in the upper mantle beneath the Chinese Tian Shan Mountains

Small-scale convection of the upper mantle beneath the Chinese Tian Shan (Tien Shan) is investigated in terms of numerical modeling. The finite element method combined with the marker-in-cell technique is used to describe the flow of the heterogeneous upper mantle. The density model is derived from the P-wave velocity structure of the crust and upper mantle along the Kuche–Kuitun profile across the Chinese Tian Shan, which is obtained using the seismic travel time tomography technique. Our computational results reveal the southward-counterclockwise and northward-clockwise upper mantle convective cells underneath the Junggar-north Tian Shan and Tarim-south Tian Shan, respectively. Our results also show the convective scale reaches to ∼500 km and the convective speed at the top of the upper mantle should not be less than 20 mm/year for a normal viscosity model. The northward extrusion of the Tarim block plays a key role in the Tian Shan mountain building since the Cenozoic period, but it nearly does not influence the upper mantle convection. The present-day tectonic deformation in the Chinese Tian Shan is related to the small-scale convection of the upper mantle.

Effects of Velocity Contrast on the Quality of Crosswell Traveltime Tomography and an Improved Method

The quality of tomographic results is not good as expected because of uneven ray coverage caused by increasing velocity contrast of neighboring geological bodies and inappropriate grid discretization in cross-well seismic first-break traveltime tomography. Tomographic results of physical and numerical models prove that when the velocity contrast is greater than 30.0%, the tomographic result will be deformed significantly, when between 30.0% and 15.0%, it will be acceptable, and when less than 15.0%, it will be very good. The multi-scale cross-well traveltime tomography is proposed on the above basis. This method is a simultaneous traveltime tomography using various schemes of grid discretization for the same model to obtain the velocity image of a study area. The multi-scale tomographic results of numerical and physical models show that this method compromises the resolution and accuracy of tomography and can significantly improve the quality of crosswell tomography. Even if the velocity contrast is more than 30.0%, the multiscale tomographic results are still good. Therefore, this method can be applied to real cross-well seismic data tomography.

Accounting for velocity anisotropy in seismic traveltime tomography: a case study from the investigation of the foundations of a Byzantine monumental building

ABSTRACTWe estimate velocity anisotropy factors from seismic traveltime tomographic data and apply a correction for anisotropy in the inversion procedure to test possible improvements on the traveltime fit and the quality of the resulting tomographic images. We applied the anisotropy correction on a traveltime data set obtained from the investigation of the foundation structure of a monumental building: a Byzantine church from the 11th century AD, in Athens, Greece.Vertical transverse isotropy is represented by one axis of symmetry and one anisotropy magnitude for the entire tomographic inversion grid. We choose the vertical direction for the symmetry axis by analysing the available data set and taking into account information on the character of the foundations of the church from the literature and past excavations. The anisotropy magnitude is determined by testing a series of values of anisotropy and examining their effect on the tomographic inversion results. The best traveltime fit and image quality are obtained with an anisotropy value (Vmax/Vmin) of 1.6, restricted to the high velocity structures in the subsurface. We believe that this anisotropy value, which is significantly higher than the usual values reported for near‐surface geological material, is related to the fabric of the church foundations, due to the shape of the individual stone blocks and the layout of the stonework.Inversion results obtained with the correction for anisotropy indicate that both the traveltime fit and the image quality are improved, providing an enhanced reconstruction of the velocity field, especially for the high‐velocity features. Based on this enhanced and more reliable reconstruction of velocity distribution, an improved image of the subsurface material character was made possible. In particular, the pattern and state of the church foundations and possible weak ground material areas were revealed more clearly. This improved subsurface knowledge may assist in a better design of restoration measures for monumental buildings such as Byzantine churches.

The crustal structure under Sanjiang and its dynamic implications: Revealed by seismic reflection/refraction profile between Zhefang and Binchuan, Yunnan

The fault belts in Sanjiang mainly include Jinshajiang-Honghe fault, Lancangjiang fault and Nujiang fault (called Sanjiang faults) in western Yunnan Province, China. By interpreting the wide-angle seismic reflection/refraction profile between Zhefang and Binchuan, which crosses Tengchong and Baoshan blocks in Dianxi (western Yunnan) tectonic zone, we reconstruct the crustal structure with seismic traveltime tomography for crustal P-wave velocity and the seismic scattering image for crustal seismic reflection structure. In this paper, we firstly present the crustal structure images of P-wave velocity and seismic reflection under the wide-angle seismic profile. These results demonstrate that, the crustal velocity structure and seismic reflection structure along the profile can be divided into 3 segments, and there is an obvious difference of crustal structure among the eastern, the western and the middle segment. Generally, crustal P-wave velocities in the Baoshan segment are 0.1―0.2 km/s slower and seismic reflection amplitudes from Moho discontinuity are stronger than the other 2 segments. In the studied area, crustal thickness is about 40 km, and shows the thickening tendency from west to east along the profile. Additionally, it can be seen that there is one strong-amplitude seismic reflection event as bright points at the depths of 8―10 km, along the segment of 80―115 km of the profile (southward of Tengchong); and seismic reflection wave-field from Moho discontinuity varies obviously along the lateral direction. Finally, we make some discussions on the crustal thickening pattern in the Sanjiang fault belt, structural environment of earthquake development and the contact rela-tionship between the Tengchong block, Banshan block and Luxi trough.

Recursive travel-time inversion: A tool for real-time seismic tomography

Abstract A new recursive inverse scheme is applied to a currently popular problem named seismic travel-time tomography, in order to enhance the efficiency and reliability in obtaining a new velocity model if a small number of new data are added to a large data set in the past. In comparison with conventional inverse schemes in seismic tomography, either least-squares or iterative types, this scheme does not require large amounts of matrix-type computations but utilizes the amount of modification in model parameters responsible for each new data set. We also introduce the computation of a collocation travel time (i.e., from a given station to every grid point) for the reference velocity model inverted by the data for all the past events, using a ray tracing scheme called the Huygens’ method (Saito, 2001), suitable to computations prior to a new event. Combining the above information already stored with the recursive inverse scheme, we can obtain a new or updated velocity model immediately after a new event takes place, because a temporal interval between two events is usually very long in a given local area. Since the model is revised at each recursive step, we perform ray tracings with the updated reference model to get more accurate ray paths and travel times than the conventional inversion schemes that use all the ray tracings for the same reference model. We first showed the validity and stability of the proposed method with synthetic data. We then applied the new approach to the P-wave travel-time data recorded in the Hidaka, south-central Hokkaido, Japan, region, and compared our result with other previous results. Our result shares the overall feature with the previous ones. In addition, a new low-velocity zone is detected in the east of the Hidaka mountains at the depth of 10 km, corresponding to the collision zone of two arcs, due to the use of the updated reference velocity model at each recursive step. We also confirmed that the order of data does not affect the final result, so that the present approach is shown as an appropriate tool for so-called real-time seismic tomography: a updated velocity model is immediately obtained at each time that a new event takes place, in order to monitor temporal variations of model parameters such as velocity structure on the real-time basis.