by GIBERT Dominique - 2 April 2006
This article presents 1-D local geo-electrical models of La Soufrière lava dome. These models were inverted with a Bayesian inversion formalism explained below.
The data used in this study are available here. If you find these data useful and use them in any publication, please cite the following paper and contact the corresponding author :
Electrical Tomography of La Soufrière of Guadeloupe Volcano: Field Experiments, 1D Inversion and Qualitative Interpretation by Florence Nicollin, Dominique Gibert, François Beauducel, Georges Boudon, and Jean-Christophe Komorowski, Earth and Planetary Science Letters, Vol. 244, 709-724, 2006.
Parameterization of the Inverse Problem Several parts of the pseudo-sections (Figure 1) acquired on La Soufrière display a 1D-like structure where the resistivity varies only in the z direction. In such cases, we invert the data under the 1D approximation to obtain vertical resistivity soundings in order to locally check for the quality of the pseudo-sections. In particular, it is possible to verify that the thicknesses of resistive layers and the values of apparent resistivity in the pseudo-sections are reasonable with respect to the inverted resistivity soundings.
Figure 1: General views of the geo-electrical pseudo-sections of La Soufrière lava dome. Regions appearing in blue correspond to low (i.e. \approx 10 \; \Omega .\textrmm) electrical conductivity and yellow is for high (i.e. \geq 5000 \; \Omega .\textrmm) resistivity. Views are for and elevation of 30 degrees. A) toward the North, B) toward the West, C) toward the South, D) toward the East.
Numerous studies have been devoted to 1D geo-electrical inversions where the unknown parameters to be determined are the thicknesses, \Delta z_n , and conductivities, \sigma _n , of the n=1,\cdots ,N layers forming the resistivity structure overlying a half-space with conductivity \sigma _0 (Gyulai 1999, Muiuane 1999, Dahlin 2001). However, in the present study we want to not only perform a 1D inversion to obtain the \Delta z’\textrms and the \sigma ’\textrms , but we also want to test the validity of the 1D approximation. Consequently, we formulate the inverse problem to be solved as follows: find the best 1D resistivity model \mathbfm with the largest data subset \mathbfd_1D compatible with the 1D approximation. The subset \mathbfd_1D counts K elements extracted from a larger data set \mathbfd corresponding to a region where the pseudo-section has a 1D appearance. In what follows, we present a non-linear inversion method where \mathbfm and \mathbfd_1D are simultaneously inverted.
The model vector is defined as
\mathbfm=\left\ \ln \sigma _0,\ln \sigma _1,\cdots ,\ln \sigma _N,\Delta z_1,\cdots ,\Delta z_N\right\ ,
with the logarithms of the conductivities taken as the unknown parameters. This parametrization is often recommended because of the large range spanned by the conductivity values. Also, the logarithmic parameterization gives an equivalent role to both the conductivity and the resistivity which are reciprocal physical quantities (Tarantola 2005). The logarithm of the conductivity is a fundamental quantity which appears in the integral equation formulation of the Poisson equation giving the perturbations of the electrical potential caused by conductivity heterogeneities (Pessel and Gibert, 2003).
The data vector, \mathbfd , consists in the logarithms of the resistance measurements R_l=\Delta V_l/I_l where \Delta V_l is the electrical potential measured between electrodes M and N , and I_l is the electrical current injected between electrodes A and B . Because of the large range spanned by the electrical potential measurements, some renormalization is necessary and, in the present study, we apply a prewithening by dividing each data by the logarithm of the squared geometrical factor deduced from the electrode arrangement. By this way, all data have almost the same magnitude and weight equivalently in the misfit estimate.
In practice, the model selection relies on the maximisation of the following likelihood function,
L\left( \mathbfm,\mathbfd_1D\right) \propto \exp\left( -\frac\left| \mathbfd_1D-\mathbfd_\mathbfm\right| _1S K^\alpha
where \left| \cdot \right| _1 stands for the L_1 norm, \mathbfd_\mathbfm is the synthetic data vector produced by the 1D conductivity model \mathbfm , K is the dimension of the data subset actually used to evaluate the fit with the model, and \alpha =3 is an exponent experimentally chosen to apply a penalty to models fitted with a relatively small number of data. By this way, the likelihood defined above allows to compare models fitted with a variable number K of data taken in the initial data set. The standard deviation, S , is initially chosen equal to unity.
Inversion by Simulated Annealing The search for the best model, \mathbfm_\textrmbest , and its corresponding data subset, \mathbfd_1D,\textrmbest , is performed with the simulated annealing algorithm which easily allows to tackle with sophisticated likelihood functions as defined above. In the following, we shall only give the technical details specifically adjusted for the present study, and the reader is referred to (Pessel and Gibert, 2003; Metropolis at al., 1953; Kirkpatrick et al., 1983; Bhanot 1988; Sambridge 2002) for general considerations concerning the Metropolis and the simulated annealing algorithms in the framework of non-linear inversion.
Let us recall that simulated annealing consists in doing Metropolis loops while applying a topological transformation to L by varying a control parameter traditionally called the temperature T>0 (Pessel and Gibert, 2003; Gibert and Virieux, 1991; Mosegaard 1995). The transformed likelihood reads,
L_T\left( \mathbfm,\mathbfd_1D\right) =L^1/T\left( \mathbfm,\mathbfd_1D\right) .
The simulated annealing loops begin with T\rightarrow \infty for which L_\infty equals a constant and finish at T=1 , i.e. at the likelihood defined above. For each temperature T , a Metropolis loop is performed to generate a sequence of models distributed according to L_T . In practice, this consists in generating a sequence of models where the next model, \left( \mathbfm,\mathbfd_1D\right) _j+1 , to be added in the sequence is obtained from the preceding one, \left( \mathbfm,\mathbfd_1D\right) _j , according to the random choice, \textrmprob\left[ \left( \mathbfm,\mathbfd_1D\right) _j+1=\left( \mathbfm,\mathbfd_1D\right) _\textrmtry\right] =\min \left[ 1,\fracL_T\left( \mathbfm,\mathbfd_1D\right) _\textrmtry
L_T\left( \mathbfm,\mathbfd_1D\right) _j\right] ,
where \left( \mathbfm,\mathbfd_1D\right) _\textrmtry is a candidate model which may eventually be included in the sequence of models. The equation above indicates that a more-likely model is always accepted and that a less-likely model is sometimes accepted. When the random assignment given by the above equation fails, i.e. when the less-likely candidate model has been rejected, the replicating transition \left( \mathbfm,\mathbfd_1D\right) _j+1=\left( \mathbfm,\mathbfd_1D\right) _j is used instead. As the Metropolis loop proceeds a sequence of models is generated such that more models fall in the most likely regions of the space model. As the temperature is lowered, these regions have their likelihood dramatically (i.e. exponentially) augmented, and provided the temperature decrease is sufficiently slow, the chain of model is gently guided toward the regions of the space model where the likelihood is maximum (Metropolis et al., 1953; Bhanot 1988).
The candidate model, \left( \mathbfm,\mathbfd_1D\right) _\textrmtry , is obtained by perturbing the current model \left( \mathbfm,\mathbfd_1D\right) _j in order to give some memory to the Metropolis chain. This point is of a particular importance in order not to make the algorithm a simple Monte Carlo search. The data subset \mathbfd_1D,\textrmtry of the candidate model is obtained by the following way: first we assign \mathbfd_1D,\textrmtry=\mathbfd_1D,j and next, in 50\% of the time, a data element is randomly chosen in the whole data set \mathbfd and is deleted from \mathbfd_1D,\textrmtry if present and is incorporated if absent. Hence, the data subset of the candidate model differs of at most one element with respect to the data subset of the last model in the Metropolis chain. The remaining parameters of the candidate model, namely the log-conductivities and the layer thicknesses, are obtained by randomly choosing a parameter, either log-conductivity or thickness, and by randomly perturbing it in a limited range.
Synthetic example We now discuss a synthetic example showing how the algorithm performs in presence of an incoherent data set whose elements were generated with two different conductivity structures. For this purpose, we invert the data set obtained by merging two synthetic data subsets of apparent resistivities computed for two different conductivity distributions (Figure 2a). In a first stage the inversion is done with no estimate of the standard deviation S which is then set to unity. The annealing schedule is performed with a temperature range 1>T>10^-4 in order to converge toward the most likely model. Figure 2b shows the likelihood of the models successively incorporated in the Metropolis chain. The resulting curve is typical of the convergence property of simulated annealing: the chain starts with randomly chosen models with a low likelihood and with strongly varying parameter values as shown in Figures 1c and 1d. As the temperature further decreases, the retained models progressively get a higher likelihood until the curve reaches a sharp step. At this stage of the cooling schedule, the model sequence suddenly gets confined in a region of high-likelihood models and the parameter values cease to strongly fluctuate (Figures 2c and 2d).
Figure 2: 1D inversion process of synthetic data: (a) the inverted pseudo-section is obtained by merging two data subsets corresponding to two models of conductivity shown on the left and on the right parts of the Figure, (b) likelihood curve of the models forming the Metropolis chain versus annealing temperature, (c) electrical conductivity versus annealing temperature, (d) top layer thickness versus annealing temperature, (e) probability distribution of the conductivities of the final model, (f) probability distribution of the top layer thickness of the final model.
The success of the global convergence of simulated annealing toward the global maximum of L is controlled by the speed of cooling which must be sufficiently low to allow an exhaustive search in the space model. Once the best model is obtained, the residuals are used to compute an estimate of S , and the last step of the simulated annealing algorithm consists in a heat-bath sequence (i.e. a Metropolis chain constructed at a constant temperature) performed at T=1 . The distribution of the models forming the heat-bath sequence may be used to construct the marginal probability curves (Gibert and Virieux, 1991; Gibert 1994; Mosegaard 1998) for each parameter as shown in Figures 1e and 1f. Depending on the random assignments made in the Metropolis algorithm, the model sequence converges toward either of the two models used to construct the synthetic data set.
Examples with real data Figure 3 shows the result of an inversion performed with a data subset extracted from the pseudo-section located on the western side of the volcano (see Figure 1). The marginal probability curves (Figures 3b and 3c) indicate a two-layer conductivity structure with a resistive layer, \rho _1\approx 3000\textrm \Omega .\textrmm and \Delta h_1\approx 35\textrm m , overlying a half-space with a resistivity \rho _0\approx 50\textrm \Omega .\textrmm . As can be observed in Figure 3d, several data have been excluded from the original data set on the right side of the pseudo-section. Indeed, these data correspond to anomalous high apparent resistivity values incompatible with the assumed 1D geometry.
Figure 3: 1D inversion results of the western part of the profile achieved around the base of the lava dome, along Chemin des Dames: (a) preselected data set, (b) probability distribution of the conductivities of the final model, (c) probability distribution of the top layer thickness of the final model, (d) data set consistent with the final model.
Another example is shown in Figure 4 for data corresponding to a fumarolic area located at the southern basis of the dome (Figure 1). In this instance, the simulated annealing inversion gives marginal probabilities such that a simple half-space model is obtained with a resistivity \rho _0\approx 20\textrm \Omega .\textrmm . Interestingly, the marginal probability curve for \Delta h_1 is flat, indicating that this parameter remains unresolved. The inverted data subset shown in Figure 4d shows that lateral data corresponding to higher apparent resistivities have been automatically eliminated.
Figure 4: 1D inversion results of a part of the profile G-G’ selected in the fumarolic area: (a) preselected data set, (b) probability distribution of the conductivities of the final model, (c) probability distribution of the top layer thickness of the final model, (d) data set consistent with the final model.
Results A total of 17 data sets corresponding to almost 1D structures recognised in the pseudo-sections of Figure 1 have been inverted according to the methodology described above. Several tests have been performed for models with more or less layers overlying the lower half-space to eventually represent a transition layer with variable resistivity. We found that models with two or more layers do not significantly improve the fit to the data. Moreover, increasing the number of layers often leads to models with low resistivity contrast and unresolved thicknesses. Consequently, according to the principle of parsimony, we decided to retain the simplest models, i.e. those with either one or no overlying layer, able to fit the data.
The results of the inversions displayed in Figure 5. All but 3 data sets give a two-layer structure with a single layer overlying a more conductive half-space. The remaining two data sets correspond to an homogeneous half-space with no variations of conductivity in the vertical direction. They correspond on the lava dome to the Jardin Lherminier area, a shallow depression with a permanent small pond formed within relatively thick phreatic ash deposits, and to active fumarolic areas characterised by a high degree of hydrothermal alteration and upflow of hydrothermal fluids located on the summit of the lava dome (Cratère Sud fumaroles) and at the South base of the lava dome near the road (Route de la Citerne fumaroles).
Figure 5: Synthesis of the 1D inversion results obtained for 17 data subsets belonging to regions where the pseudo-sections appear laterally invariant. The great rectangle represents the resistivity of the bottom half space, the line represents the resistivity of the top layer and its thickness is proportional to the top layer thickness. The values of altitude correspond to the depth of the top of the conductive half space.
Conclusions Figure 5 reveals that the 17 geo-electrical soundings constitute a simple global conductivity structure with a low-resistivity (\rho _0\approx 50\textrm \Omega .\textrmm) basal layer on the eastern and southern edges of the lava dome, and a medium-resistivity (\rho _0\approx 250\textrm \Omega .\textrmm) basal layer on the northern edge and on the lava dome itself. In most cases, the basal layer is overlained by a resistive (\rho _1\approx 2000\textrm \Omega .\textrmm) layer with a thickness \Delta h_1 comprised between 10 and 30 meters. Noticeable exceptions are for the August 30th Fracture and the Carbet areas where the overlying layer has a significantly higher resistivity (\rho _1\approx 4000\textrm \Omega .\textrmm) and corresponds to rock avalanche blocks. Another exception is for the Col de l’échelle area where the overlying layer has a lower resistivity (\rho _1\approx 400\textrm \Omega .\textrmm) and corresponds to fumarolic zones that ceased to be active in 1979 and 1984 (Boudon 1989,Komorowski 2005,Komorowski 2001).
Acknowledgments This work is dedicated to the memory of our Friend and Colleague Alberto Tarchini who participated to most of our field experiments and always shared his generous enthusiasm. This study benefited from invaluable help from the whole staff of the Observatoire Volcanologique et Sismologique de Guadeloupe: Christian Anténor-Habazac, Sara Bazin, Véronique Daniel, Bertrand Figaro, Gilbert Hammouya, Christian Lambert, Didier Mallarino and Laurent Mercier. D.G. and F.N. would like to thank Gérard Werther for numerous discussions concerning historical records about La Soufrière. Financial support was provided through the CNRS/INSU Antilles Program and from IPGP .
Cited references D. Gibert, M. Pessel, Identification of sources of potential fields with the continuous wavelet transform: Application to self-potential profiles, Geophys. Res. Lett. 28 (2001) 1863-1866. K. Aizawa, R. Yoshimura, N. Oshiman, K. Yamazaki, T. Uto, Y. Ogawa, S.B. Tank, W. Kanda, S. Sakanaka, Y. Furukawa, T. Hashimoto, M. Uyeshima, T. Ogawa, I. Shiozaki, A.W. Hurst, Hydrothermal system beneath Mt. Fuji volcano inferred from magnetotellurics and electric self-potential, Earth and Planet. Sci. Lett. in the press (2005). A. Gyulai, T. Ormos, A new procedure for the interpretation of VES data: 1.5D simultaneous inversion method, J. Applied Geophys. 41 (1999) 1-17. E.A. Muiuane, L.B. Pedersen, Automatic 1D interpretation of DC resistivity sounding data, J. Applied Geophys. 42 (1999) 35-45. T. Dahlin, The development of DC resistivity imaging techniques, Computers Geosciences 27 (2001) 1019-1029. A. Tarantola, Elements for Physics: Quantities, Qualities, and Intrisic Theories, Springer, Berlin, (2006) 263 pp. M. Pessel, D. Gibert, Multiscale electrical impedance tomography, J. Geophys. Res. 108B1 (2003) 2054, doi:10.1029/2001JB000233. N. Metropolis, A. Rosenbluth, N. Rosenbluth, A. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087-1092. S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671-680. G. Bhanot, The Metropolis algorithm, Rep. Prog. Phys. 51 (1988) 429-457. M. Sambridge, K. Mosegaard, Monte Carlo methods in geophysical inverse problems, Rev. Geophys. 40(3) (2002) 1009, doi:10.1029/2000RG00089. D. Gibert, J. Virieux, Electromagnetic imaging and simulated annealing, J. Geophys. Res. 96 (1991) 8057-8067. K. Mosegaard, A. Tarantola, Monte-Carlo sampling of solutions of inverse problems, J. Geophys. Res. 100 (1995) 12431-12447. D. Gibert, B. Tournerie, J. Virieux, High-resolution electromagnetic imaging of the conductive Earth interior, Inv. Problems 10 (1994) 341-351. K. Mosegaard, Resolution analysis of general inverse problems through inverse Monte Carlo sampling, Inv. Problems 14 (1998) 405-426.