Skip to main content

Markov random field modeling for mapping geofluid distributions from seismic velocity structures

Abstract

We applied the Markov random field model, which is a kind of a Bayesian probabilistic method, to the spatial inversion of the porosity and pore shape in rocks from an observed seismic structure. Gaussian Markov chains were used to incorporate the spatial continuity of the porosity and the aspect ratio of the pore shape. Synthetic inversion tests were able to show the effectiveness and validity of the proposed model by appropriately reducing the statistical noise from the observations. The proposed model was also applied to natural data sets of the seismic velocity structures in the mantle wedge beneath northeastern Japan, under the assumptions that the fluid was melted and the temperature and petrologic structures were uniformly distributed. The result shows a significant difference between the volcanic front and the forearc regions, at a depth of 40 km. Although the parameters and material properties will need to be determined more precisely, the Markov random field model presented here can serve as a basic inversion framework for mapping geofluids.

Correspondence/findings

Introduction

In order to understand the various dynamic processes in the earth, it is important to understand the distribution of geofluids. Recent developments in the technology for geophysical observations, such as seismic tomography and geomagnetic methods, provide detailed images of the earth’s interior (Nakajima et al. 2001; Ogawa et al. 2001; Takahashi et al. 2009). Additionally, there has been increased understanding of the constitutive relationships between the physical variables, such as lithology, the porosity of rocks, and the observational data, such as seismic velocity and resistivity (Glover et al. 2000; Takei 2002). Against this background, a pioneering study by Nakajima et al. (2005) used the constitutive function proposed by Takei (2002) to evaluate the effective aspect ratio and the volume fraction (porosity) of the fluid-filled pores in the observed low-velocity anomalies in the mantle wedge beneath northeastern Japan.

Recently, several studies have attempted to make a quantitative and detailed map of the spatial distribution of geofluids (Hoshide and Nakamura 2013; Iwamori et al. 2011). However, this remains difficult, because there is still much uncertainty in the available data and assumptions. In order to overcome the difficulties arising from this noise and uncertainty, a statistical and probabilistic analysis of the geophysical data is essential.

The main purpose of the present study is to construct an inversion framework that can be used to estimate precisely the distributions of various physical properties from observed spatial data sets; we do this by developing the Markov random field (MRF) model, which is a kind of a Bayesian statistical model. The Bayesian approach enables us to incorporate a forward model and prior information into a data-driven inversion analysis.

The MRF model uses Markov chains to describe the properties of an image, and it is often used in the field of information science for image restoration and pattern recognition (Geman and Geman 1984; Li 2009; Tanaka 2002). In the MRF model, the spatial variations in physical properties are assumed to be generally smaller than the noise in the data and the analytical uncertainty. If this assumption is valid, then by using the Bayesian approach, the MRF model appropriately filters out the high-frequency noise, and we can obtain the accurate spatial distributions of the physical properties. Recent papers in the natural sciences have applied the MRF model to inversion problems for various observational data sets (Kuwatani et al. 2012; Watanabe et al. 2009).

Here, we develop a Gaussian MRF model to reconstruct the spatial distribution of geofluids from the seismic velocity structure. On the basis of the Bayesian framework, the process for generating the velocity structure and the spatial continuity of the distribution of geofluids are introduced into the stochastic inversion analysis in accordance with the law of causality. In order to deal with the nonlinear relationship between the target physical variables and the observed data, a Markov chain Monte Carlo (MCMC) algorithm was incorporated into the MRF model (Metropolis et al. 1953). An application of the method to synthetic data showed that the spatial distributions of porosity and the aspect ratio could be reliably estimated, and this supports the effectiveness of the MRF model. We also applied the model to the velocity structure of the mantle wedge beneath northeastern Japan, which was obtained by 3-D tomography (Matsubara et al. 2008), under the simple assumption that variables other than the porosity and aspect ratio were known and spatially uniform. Finally, we will discuss the validity of our assumptions, the effectiveness and applicability of the MRF model, and the geophysical implications. Although many parameters and material properties remain to be determined more precisely, the proposed framework will be very effective for determining the distribution of geofluids.

Method

The seismic wave velocities (VP and VS) of solid-fluid composite media are generally expressed as functions of the intrinsic elastic parameters of the solid framework and pore-filling fluid, fluid volume fraction, and pore geometry (Mavko 1980; Takei 2002). Recently, (Takei 2002) proposed a unified formulation of VP and VS as a function of the effective aspect ratio (α) and the fluid volume fraction, that is, the porosity (ϕ), of the fluid-filling pores; this formulation can be applied to a wide variety of pore shapes [see Additional file 1]. Assuming that the type and compositions of rock and fluid are fixed and that the thermal and pressure effects are negligible, these functions can be simply rewritten as follows:

V P = f P ( ϕ , α ) , V S = f S ( ϕ , α ) .
(1)

As shown in Figure 1, seismic velocities show a monotonic decrease with increasing porosity ϕ, and the slope is controlled by the aspect ratio α.

Figure 1
figure 1

Dependency of seismic wave velocities on porosity and effective aspect ratio, calculated using Takei (2002). (a) P-wave velocity and (b) S-wave velocity. The dry rock and fluid phase are assumed to be peridotite ( V P 0 =7.9 km/s, V S 0 =4.55 km/s) and melt.

However, geophysical observations always contain some measurement errors and uncertainty. Thus, the obtained P-wave velocity V P i for each grid cell i should be written as

V P i = f P ( ϕ i , α i )+ ε P i ,
(2)

where ε P i is the observational noise for V P i for each spatial grid cell (i). If we assume a Gaussian noise with zero mean, Equation 2 can be rewritten in terms of the conditional probability as

p( V P i | ϕ i , α i )= 1 2 π σ P 2 exp - V P i - f P ( ϕ i , α i ) 2 2 σ P 2 ,
(3)

where p( V P i | ϕ i , α i ) is the probability that V P i is generated, given ϕi and αi, and σ P 2 is the variance of the noise in the observed seismic velocity VP. Equations 2 and 3 can also be written for the S-wave seismic velocity VS. Because we will assume that random errors are independent, the probabilities of generating VP and VS are independent between grid cells. By multiplying the probabilities of VP and VS for all grid cells, we can obtain the total joint probability for the observed velocities VP and VS for all grid cells as

p( V P , V S |ϕ,α)= i = 1 N p( V P i | ϕ i , α i )·p( V S i | ϕ i , α i ),
(4)

where N is the total number of grid cells measured, and V P , V S , ϕ, and α indicate the respective set of variables VP, VS, ϕ, and α for the observed grid cells i=1,…N.

On the other hand, Bayes’ theorem can be written as follows:

p(ϕ,α| V P , V S )= p ( V P , V S | ϕ , α ) · p ( ϕ , α ) p ( V P , V S ) .
(5)

The left-hand side of the equation is the posterior probability p(ϕ,α|V P ,V S ), where the probability of ϕ and α is based on the observed seismic velocities, V P and V S . The set of ϕ and α that maximizes the posterior probability is considered to be the most probable spatial distribution of porosity ϕ and aspect ratio α for the available observed data set and the available prior information about the distribution of the fluid. The numerator of the right-hand side of the equation is the product of the likelihood function p(V P ,V S |ϕ,α) and the prior probability p(ϕ,α). The likelihood function p(V P ,V S |ϕ,α) is the probability of generating the observed data V P and V S , given that values ϕ and α are true. The prior probability p(ϕ,α) is the probability of the physical variables ϕ and α before any additional observations are made. We assume that the physical variables are continuous, and so the physical quantities are similar at neighboring spaces and times. The MRF model adopts the Gaussian Markov chain model as the prior probability, as follows:

p(ϕ)= 1 Z ϕ exp - 1 2 σ ϕ 2 i j ϕ i - ϕ j 2 ,
(6)

where i j is the summation of all pairs of neighboring grid cells, σ ϕ 2 is the variance of the change in ϕ between two adjacent grid cells, and Z ϕ ( σ ϕ 2 ) is a normalization coefficient. The prior probability for α can be also written as Equation 6. The denominator of the right-hand side of the equation p(V P ,V S ) is invariant for changes in ϕ and α, so this is negligible for our analysis.

Here, we define the evaluation function as the negative logarithm of the posterior possibility, - lnp(ϕ,α|V P ,V S ). By substituting the likelihood function Equation 4 and the prior probability Equation 6 into Bayes’ theorem Equation 5, the evaluation function can be expressed as

E ( ϕ , α ; θ , V P , V S ) = 1 2 σ P 2 i = 1 N V P i - f P ( ϕ i , α i ) 2 + 1 2 σ S 2 i = 1 N V S i - f S ( ϕ i , α i ) 2 + 1 2 σ ϕ 2 i j ϕ i - ϕ j 2 + 1 2 σ α 2 i j α i - α j 2 + N 2 ( ln σ P 2 + ln σ S 2 ) + ln Z ϕ + ln Z α + C ,
(7)

where θ indicates the set of parameters { σ P 2 , σ S 2 , σ ϕ 2 , σ α 2 }, and C is a constant that is independent from ϕ, α, and θ. Due to the monotonicity of the logarithm function, the minimization of the evaluation function E(ϕ,α;θ) is equivalent to the maximization of the posterior probability p(ϕ,α|V P ,V S ).

In the evaluation function, the first and second terms indicate the reproducibility of the observation, respectively, whereas the third and forth terms indicate the spatial continuity of the porosity and the aspect ratio, respectively. Minimization of E(ϕ,α;θ) satisfies the requirements of both the reproducibility of the observed data and the spatial continuity of the physical variables.

The parameters θ fully control the relative importance of the reproducing the observational data to honoring the continuity of the physical properties, and so these parameters are often referred to as hyperparameters. The most probable set of hyperparameters is obtained by minimizing the free energy, which is defined as the negative logarithm of the posterior possibility p(θ|V P ,V S ). Using Bayes’ theorem and marginalizing the likelihood function, the free energy can be expressed as

F ( θ ) - ln p ( θ | V P , V S ) = - ln - - exp - E ( ϕ , α ; θ ) d ϕ d α + C ,
(8)

where we assume that the prior probability p(θ) is uniformly distributed, and C is a constant that is independent from θ. In this study, the free energy F(θ) was minimized by the steepest descent method, using the MCMC method [see Additional file 1]. A maximum a posteriori (MAP) solution set of ϕ and α can also be obtained from numerous candidates which are generated by the MCMC calculations.

Synthetic inversion test

The synthetic inversion test was conducted to investigate how well the proposed method could reconstruct the target physical quantities from a noisy data set. The target spatial distributions of ϕ and α were assumed to have Gaussian inhomogeneities (Figure 2a). In a Gaussian inhomogeneity, the difference between the values for a physical property between two adjacent grid cells obeys a Gaussian distribution and corresponds to the prior probability of a Gaussian Markov chain (Equation 6). In order to simulate actual observed data, the velocity structure, V P and V S , was generated by substituting the assumed 2-D spatial distributions of ϕ and α into the constitutive function Equation 1 and then adding noise (Figure 2b). For the constitutive function Equation 1, we used the same physical properties as were used for the melted peridotite system: 1,050°C and a depth of 40 km (Nakajima et al. 2005).

Figure 2
figure 2

Synthetic data used in the inversion test. (a) Synthetic distributions of porosity ϕ and aspect ratio α, which are to be estimated. They were generated by 50×50 random walks; the variances were set to ( σ ϕ 2 , σ α 2 )=(0.00 4 2 ,0.01 5 2 ). (b) Observational data of V P and V S . They were obtained using Equation 1 with Gaussian noise; the variances were set to ( σ P 2 , σ S 2 )=(0. 1 2 ,0.0 5 2 ).

Figure 3 shows the changes in the hyperparameters θ during the iterations of the steepest descent method as it minimized the free energy F(θ) by using the MCMC method. We can see that each hyperparameter converges to the true value as the number of trials increases; this indicates that the estimations of the hyperparameters are successful. Figure 4 shows the root-mean-square (RMS) errors (the square root of the average of the squared residuals between the true and estimated values) for the values of ϕ and α that minimize the evaluation function E(ϕ,α;θ) for each iteration of the steepest descent method. They decrease and converge to a low value as the number of trials increases, and this implies that the estimation accuracy increases as the hyperparameter estimation proceeds.

Figure 3
figure 3

Estimated behavior of the hyperparameters during the method of steepest descent. (a) Variances of continuity of porosity (ϕ) and aspect ratio (α), and (b) variances of noise in the observational data, VP and VS. The empty circles indicate the true values of the variances.

Figure 4
figure 4

Behavior of the root-mean-square errors in the porosity ( ϕ ) and the aspect ratio ( α ). The trial number is the iteration number for the method of steepest descent.

Figure 5a shows the spatial distributions estimated by the MRF model for ϕ and α from the above synthetic velocity structure model. The distributions of ϕ and α were calculated by numerically solving Equation 1 from the observed V P and V S ; we shall refer to this as the deterministic method. For comparison, these are also shown in Figure 5b. The distributions of ϕ and α calculated by the deterministic method are directly affected by the noise of the observed seismic velocities. They are too jagged for the true distributions to be determined. Although the deterministically estimated distributions of ϕ and α are rough and contain large errors, the MAP solutions that were obtained by the MRF model are much smoother and more accurate. For the deterministic method, the RMS error for α is approximately 0.31, whereas the RMS error for ϕ is approximately 0.07. With the MRF model, the RMS errors of ϕ and α were reduced to 0.0098 and 0.038, respectively. Although the deviations of α from the true profile are still large, the estimated profiles are much more accurate than those obtained by the deterministic method.

Figure 5
figure 5

Estimated fluid distribution (porosity ϕ and aspect ratio α ). (a) MRF model (this study) and (b) deterministic method.

In addition to the synthetic distributions of ϕ and α, which have Gaussian inhomogeneities, we also checked the effectiveness of the proposed method using a power-law inhomogeneity, since this is considered to be the typical distribution of actual inhomogeneities in the earth (Sato et al. 2012). The details are in Additional file 1. Although there is difference between the prior probabilities, which were assumed to be Gaussian, and the actual power-law inhomogeneities, the estimated variances of continuity, σ ϕ 2 and σ α 2 , were approximately the same as the true values obtained by estimating them from the hyperparameters. We can also estimate the spatial distributions of ϕ and α, which supports the effectiveness of the proposed method for actual non-Gaussian distributions in natural systems. Although further investigation is needed, due to the versatility of the Gaussian distributions, the proposed method is approximately valid for natural continuous distributions.

Application to the mantle wedge beneath northeastern Japan

We applied the above MRF model to the velocity structure beneath northeastern Japan, as constructed by Matsubara et al. (2008) using seismic tomography. In Matsubara et al. (2008), the cells were constructed with a 0.1° grid spacing in the horizontal direction and a 10-km grid spacing in the vertical direction. In this study, the 2-D distributions of the fluid were imaged from the horizontal VP and VS images at a depth of 40 km (Figure 6). As in Nakajima et al. (2005), we assumed that the host rock was peridotite; for the reference velocities, defined as the velocities of the dry host rock, we used V P 0 =7.90 and V S 0 =4.55 km/s. We also assumed that the fluid was melted throughout. We used the same melt parameters as those used in Nakajima et al. (2005).

Figure 6
figure 6

Observed tomographic data ( Matsubara et al. 2008) at the depth of 40 km used in this study. P-wave and S-wave velocities (VP (left) and VS (right)) and their ratio (VP/VS, bottom) are shown. Closed triangles indicate the Quaternary volcanoes.

We show the results for the inversion beneath northeastern Japan. Figure 7 shows the estimated hyperparameters versus the iteration count during the method of steepest descent. From Figure 7, we can confirm the stable convergence of all the hyperparameters. We also checked the stability of convergence for different initial values of the hyperparameters and verified that there were no local minima that could trap the inversion. The variances of the continuity of the porosity σ ϕ 2 and the aspect ratio σ α 2 converged to 8.3×10-7 and 2.7×10-5, respectively. On the other hand, the variances of the observational noises VP and VS converge to 0.061 and 0.0038, respectively. The standard deviations σ ϕ and σ α equaled the RMS errors. The calculated RMS residuals of VP and VS were 0.25 and 0.062, respectively.

Figure 7
figure 7

Estimated behavior of the hyperparameters during the method of steepest descent, at the depth of 40 km. (a) Variance of continuity of porosity σ ϕ 2 and aspect ratio σ α 2 ; (b) variance of the observational noise in VP σ P 2 and VS σ S 2 .

Figure 8a shows the spatial distributions of ϕ and α estimated by the MRF model. The most probable estimate is the set of ϕ and α that minimizes the evaluation function E(ϕ,α;θ), as defined by Equation 7, i.e., the set that maximizes p(ϕ,α|V P ,V S ). For comparison, the distributions of ϕ and α, which were estimated by the deterministic method, are shown in Figure 8b. The distributions of ϕ and α calculated by the deterministic method are very jagged, which directly reflects the noise in the observed seismic velocities. Additionally, many of the observed V P and V S deviate from the range of possible V P and V S , as generated by Equation 1; this is as shown in Figure 8b.

Figure 8
figure 8

Estimated distributions of ϕ (left) and α (right). (a) MRF model and (b) deterministic method. The gray pixel in the estimate of the deterministic method indicates that the residual of the calculated and observed VP and VS was larger than 1% of the observed VP and VS.

The ϕ values estimated by the MRF model range from 0 to 0.01. By comparison with the original mappings of VP and VS (Figure 6), the regions of large ϕ values correspond to small VP and VS, reflecting that both VP and VS decrease monotonically with increasing ϕ. In particular, the ϕ values are relatively large (>0.002) beneath the Quaternary volcanoes and Hokkaido. On the forearc side, the value of ϕ is generally low, ranging from 0 to 0.002. However, several anomalous, large values of ϕ are found in the east, off Fukushima and between the Honshu and Hokkaido islands.

The α values are 0.01 to 0.03 on the back-arc side and 0.001 to 0.01 on the forearc side. The regions of small α are roughly consistent with the regions that have a high VP/VS ratio. Beneath the Quaternary volcanoes, the α value is generally high. In particular, near Chubu and west Hokkaido, the value may be as large as 0.1. On the forearc side, the α value is generally low, and small values of α, about 0.001, are detected in the Hidaka region and westward.

Discussion

In terms of reducing high-frequency noise, the role and efficiency of the MRF model appear to be similar to those of a smoothing filter applied to the observational data. In a smoothing method such as a moving average, too large of a filter will cause excessive smoothing and blur the details of the image. In actual analyses for natural systems, however, the true distribution and magnitude of the noise are unknown, so we cannot make a prior determination of the appropriate filter size. Thus, it is difficult to use an ordinary filtering (averaging) method to determine the precise physical values from noisy observations. However, even without a priori information about the magnitude of the noise, the MRF model can determine the variance of the noise from the data. This is the most significant advantage of the MRF model, that it enables us to analyze the data objectively and quantitatively.

When applied to the actual data, the deterministic method, which uses the data and the inverse function to obtain an analytic solution, results in a very jagged and incomplete estimate. When observational data is converted to physical parameters, the results are sometimes beyond the scope of the model, and thus no solution can be derived. This is caused by a combination of observational noise and uncertainty, which highlights the importance of using a statistical or probabilistic analysis. The proposed method was able to image the continuous distribution of fluid because it did not take a deterministic approach but a probabilistic approach, and it was thus able to avoid perturbations due to noise.

This study used the seismic velocities estimated by tomography, V P and V S , and realistic observations were simulated by adding uncorrelated Gaussian noise with zero mean to each of the grid cells. At present, the proposed model cannot deal with the cross-correlated errors of seismic velocities that are derived from a tomographic inversion. For a more accurate estimation of the distribution of geofluids, a ray-path matrix, which relates the travel time to the inverse of the seismic velocity, can be incorporated into the probabilities of the generated observations, Equation 3. In our current research, due to the flexibility of the MRF model, we were able to successfully apply it to a seismic tomographic inversion by using the ray path matrix to generate the probabilities and continuity of the rate to obtain the prior probability. The inversion of physical properties directly from the observed travel time is an important issue that needs to be addressed in future work.

The calculated spatial variations in the porosity ϕ and the aspect ratio α of the geofluids show a significant difference between the forearc regions and the volcanic front, at a depth of 40 km. On the forearc side, the values are generally low, with ϕ ranging from 0 to 0.2 vol.%, and α ranging from 0.001 to 0.01; this indicates that the fluid is not intergranular but is between thin cracks. There is little melting in this region, and even if melt exists, it is not textually equilibrated to the surrounding rocks. The small amount of melt is consistent with other geophysical observations which indicate weak inhomogeneities and weak attenuation on the forearc side (Takahashi et al. 2009; Umino and Hasegawa 1984; Yoshimoto et al. 2006).

Beneath the Quaternary volcanoes, on the other hand, the large amount of geofluid (> 0.2 vol.%) indicates partial melting of the rock. The α values are generally high (> 0.02), so the melt is considered to fill oblate spheroid cracks or dikes. In particular, in the region of Chubu and west Hokkaido, the value is up to 0.1, which indicates that the pore geometry is near equidimensional, and the fluid is distributed in the spaces between the grains (Waff and Bulau 1979; Takei 2002). The estimated amount and shape of geofluid are considered to be closely related to the magmatic process of the Quaternary volcanoes (Tamura et al. 2002; Nakajima et al. 2005).

At the same depth beneath the volcanic front, Nakajima et al. (2005) estimated the porosity at 1 to 2 vol.% with an aspect ratio of 0.02 to 0.04. Although the values of the aspect ratios are consistent with each other, the porosities differ. It is not possible to make a simple comparison between this study and that of (Nakajima et al. 2005), because different seismic velocity structures were used in their analyses (Nakajima et al. 2005; Matsubara et al. 2008). There are also several sources of uncertainty in natural systems, which affect the values of the parameters used in the analysis; this is discussed below. Further studies are necessary to evaluate the validity of the obtained distributions of geofluids.

In order to match the number of unknown parameters to the number of observable parameters, we have assumed that the parameters other than porosity and aspect ratio of the geofluids are known and uniformly distributed. In the actual mantle wedge zone, however, the spatial variations of temperature and composition overlap with geofluid distributions. For more realistic imaging of the distribution of geofluids, it is necessary to introduce a priori information and to introduce other models, such as those for thermal or petrological structure, into the analysis (Iwamori et al. 2011). However, many of these models are still poorly constrained, and thus, in order to obtain reliable distributions of geofluids, some of the parameters should be probabilistic variables. The MRF model may have the potential to overcome these difficulties due to its Bayesian approach and flexible formalism. By allowing us to add terms to the evaluation function, it allows us to incorporate other geophysical observational data sets and various types of prior information as probabilistic constraints. Although additional theoretical improvements are needed for individual problems, the MRF model presented here can serve as a basic inversion framework for the mapping of geofluids.

References

  • Geman S, Geman D: Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE T Pattern Anal 1984, 6(6):721–741.

    Article  Google Scholar 

  • Glover PWJ, Hole MJ, Pous J: A modified Archie’s law for two conducting phases. Earth Planet Sci Lett 2000, 180(3–4):369–383.

    Article  Google Scholar 

  • Hoshide T, Nakamura M: Development of “Geofluid Map” in the Naruko district, NE Japan. 2013. Paper presented at the Japan Geoscience Union meeting, Chiba, Japan, 19–24 May 2013 (SCG63-P10)

    Google Scholar 

  • Iwamori H, Watanabe T, Nakamura M, Ichiki M, Nakajima J, Ogawa Y, Okada T, Matsuzawa T: An integrated model for mapping geofluids. 2011. Paper presented at the American Geophysical Union fall meeting, San Francisco, 5–9 Dec 2011 (V41D-2545)

    Google Scholar 

  • Kuwatani T, Nagata K, Okada M, Toriumi M: Precise estimation of pressure-temperature paths from zoned minerals using Markov random field modeling: theory and synthetic inversion. Contrib Mineral Petrol 2012, 163(3):547–562. 10.1007/s00410-011-0687-3

    Article  Google Scholar 

  • Li SZ: Markov random field modeling in image analysis. Heidelberg: Springer; 2009.

    Google Scholar 

  • Matsubara M, Obara K, Kasahara K: Three-dimensional P- and S-wave velocity structures beneath the Japan islands obtained by high-density seismic stations by seismic tomography. Tectonophysics 2008, 454(1–4):86–103. 10.1016/j.tecto.2008.04.016

    Article  Google Scholar 

  • Mavko GM: Velocity and attenuation in partially molten rocks. J Geophys Res 1980, 85(10):5173–5189.

    Article  Google Scholar 

  • Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E: Equation of state calculations by fast computing machines. J Chem Phys 1953, 21(6):1087–1092. 10.1063/1.1699114

    Article  Google Scholar 

  • Nakajima J, Matsuzawa T, Hasegawa A, Zhao DP: Three-dimensional structure of V p , V s , and V p / V s beneath northeastern Japan: implications for arc magmatism and fluids. J Geophys Res Solid Earth 2001, 106(B10):21843–21857. 10.1029/2000JB000008

    Article  Google Scholar 

  • Nakajima J, Takei Y, Hasegawa A: Quantitative analysis of the inclined low-velocity zone in the mantle wedge of northeastern Japan: a systematic change of melt-filled pore shapes with depth and its implications for melt migration. Earth Planet Sci Lett 2005, 234(1–2):59–70. 10.1016/j.epsl.2005.02.033

    Article  Google Scholar 

  • Ogawa Y, Mishina M, Goto T, Satoh H, Oshiman N, Kasaya T, Takahashi Y, Nishitani T, Sakanaka S, Uyeshima M, Honkura Y, Matsushima M: Magnetotelluric imaging of fluids in intraplate earthquake zones, NE Japan back arc. Geophys Res Lett 2001, 28(19):3741–3744. 10.1029/2001GL013269

    Article  Google Scholar 

  • Sato H, Fehler M, Maeda T: Seismic wave propagation and scattering in the heterogeneous Earth. New York: Springer; 2012.

    Book  Google Scholar 

  • Takahashi T, Sato H, Nishimura T, Obara K: Tomographic inversion of the peak delay times to reveal random velocity fluctuations in the lithosphere: method and application to northeastern Japan. Geophys J Int 2009, 178: 1437–1455. doi:10.1111/j.1365–246X.2009.04227.x

    Article  Google Scholar 

  • Tamura Y, Tatsumi Y, Zhao D, Kido Y, Shukuno H: Hot fingers in the mantle wedge: new insights into magma genesis in subduction zones. Earth Planet Sci Lett 2002, 197: 105–116. 10.1016/S0012-821X(02)00465-X

    Article  Google Scholar 

  • Takei Y: Effect of pore geometry on V P / V S : from equilibrium geometry to crack. J Geophys Res Solid Earth 2002, 107(B2):2043–2054. doi:10.1029/2001JB000522 doi:10.1029/2001JB000522

    Article  Google Scholar 

  • Tanaka K: Statistical-mechanical approach to image processing. J Phys A Math Gen 2002, 35(37):R81. doi:10.1088/0305–4470/35/37/201

    Article  Google Scholar 

  • Umino N, Hasegawa A: Three-dimensional Qs structure in the northeastern Japan arc. J Seismol Soc Japan Ser 1984, 37(2):217–228. (in Japanese with English abstract) (in Japanese with English abstract)

    Google Scholar 

  • Waff HS, Bulau JR: Equilibrium fluid distribution in an ultramafic partial melt under hydrostatic stress conditions. J Geophys Res 1979, 84: 6109–6114. 10.1029/JB084iB11p06109

    Article  Google Scholar 

  • Watanabe K, Tanaka H, Miura K, Okada M: Transfer matrix method for instantaneous spike rate estimation. IEICE Trans Inform Sys E92D 2009, 7: 1362–1368.

    Article  Google Scholar 

  • Yoshimoto K, Wegler U, Korn A: A volcanic front as a boundary of seismic-attenuation structures in northeastern Honshu, Japan. Bull Seismol Soc Am 2006, 96: 637–646. doi:10.1785/0120050085

    Article  Google Scholar 

Download references

Acknowledgements

We would like to thank the editor and two anonymous reviewers for their valuable comments. Discussions with H. Iwamori and A. Okamoto were thoughtful and constructive. We used the seismic velocity structure of (Matsubara et al. 2008), provided by the National Research Institute for Earth Science and Disaster Prevention (NIED). This study was financially supported by the research project ‘Evaluation and disaster prevention research for the coming Tokai, Tonankai, and Nankai earthquakes’ from the MEXT of Japan. This work was also supported by JSPS KAKENHI Grant Numbers 25120005, 25120009, and 25280090.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Tatsu Kuwatani.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

TK, MO, and MT designed the research. TK and KN performed the research. KN and MO contributed the new algorithm. TK wrote the paper. All authors read and approved the final manuscript.

Electronic supplementary material

40623_2013_5_MOESM1_ESM.pdf

Additional file 1:Markov random field modeling for mapping geofluid distributions from seismic velocity structures by T. Kuwatani et al.(PDF 191 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Kuwatani, T., Nagata, K., Okada, M. et al. Markov random field modeling for mapping geofluid distributions from seismic velocity structures. Earth Planet Sp 66, 5 (2014). https://doi.org/10.1186/1880-5981-66-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1880-5981-66-5

Keywords