Time-resolved tomography using acoustic emissions in the laboratory, and application to sandstone compaction

S U M M A R Y Acoustic emission (AE) and active ultrasonic wave velocity monitoring are often performed during laboratory rock deformation experiments, but are typically processed separately to yield homogenized wave velocity measurements and approximate source locations. Here, I present a numerical method and its implementation in a free software to perform a joint inversion of AE locations together with the 3-D, anisotropic P-wave structure of laboratory samples. The data used are the P-wave first arrivals obtained from AEs and active ultrasonic measurements. The model parameters are the source locations and the P-wave velocity and anisotropy parameter (assuming transverse isotropy) at discrete points in the material. The forward problem is solved using the fast marching method, and the inverse problem is solved by the quasi-Newton method. The algorithms are implemented within an integrated free software package called FaATSO (Fast Marching Acoustic Emission Tomography using Standard Optimisation). The code is employed to study the formation of compaction bands in a porous sandstone. During deformation, a front of AEs progresses from one end of the sample, associated with the formation of a sequence of horizontal compaction bands. Behind the active front, only sparse AEs are observed, but the tomography reveals that the P-wave velocity has dropped by up to 15 per cent, with an increase in anisotropy of up to 20 per cent. Compaction bands in sandstones are therefore shown to produce sharp changes in seismic properties. This result highlights the potential of the methodology to image temporal variations of elastic properties in complex geomaterials, including the dramatic, localized changes associated with microcracking and damage generation.


I N T RO D U C T I O N
In situ imaging of samples during laboratory rock deformation experiments is crucial for our understanding of the physics of deformation and strain localization in rocks. However, the elevated pressures and temperatures required to perform relevant experiments generally preclude direct observation techniques. Since the 1960s, in situ acoustic emission (AE) monitoring techniques have been developed to detect microfracture propagation and extract event locations, amplitudes and focal mechanisms (e.g. Mogi 1968;Scholz 1968;Nishizawa et al. 1984;Lockner 1993); such observations have revealed unprecedented details about the dynamics of fault growth (Lockner et al. 1992) and strain localization in crustal rocks (e.g. Fortin et al. 2009). Passive AE monitoring is only able to record the microcracks that propagate dynamically and generate elastic waves within the frequency and amplitude range of the recording devices. Therefore, passive AE monitoring can be usefully complemented by active measurements of P-and S-wave velocities, which provide unique information on the development of seismic anisotropy (e.g. Gupta 1973;Lockner & Byerlee 1977) and on the evolution of microfracture density and orientation (e.g. Sayers & Kachanov 1995;Schubnel et al. 2003) during rock deformation experiments.
Despite the great advances in in situ imaging made possible by AE and ultrasonic monitoring, both methods rely on the key simplifying assumption that the elastic wave velocities remain homogeneous throughout the material. Indeed, AE locations are determined by assuming a homogeneous, possibly anisotropic velocity structure and the active measurements of P-and S-wave velocities are made between pairs of sensors assuming straight ray paths, effectively corresponding to path-averaged velocities. However, it is well known that strong structural heterogeneities often develop in rock samples during deformation experiments, especially in the brittle regime where strain and microcrack damage localize along a fault plane. We expect these strain and damage heterogeneities to be associated with sharp local changes in elastic wave velocities and anisotropy, by up to several tens of percent (e.g. Schubnel 2178N. Brantut et al. 2003Fortin et al. 2006;Brantut et al. 2014). The development of a complex wave velocity structure in the material has two consequences. First, if not properly accounted for, it severely limits the accuracy of AE locations. Secondly, and perhaps most importantly, the wave velocity structure becomes a key signature of the deformation processes operating in the material. It is therefore tempting to aim to determine a tomographic picture of the velocity structure, using AE and/or ultrasonic data.
In order to achieve this, the most natural starting point is to use active measurements with known source positions. Using this type of data, Yukutake (1989), Falls et al. (1992) and Jansen et al. (1993) were able to retrieve 2-D tomographic sections of rock samples subjected to shear, thermal and hydraulic stresses, respectively. In all cases, the array of transducers was quite extensive, the rock samples were quite large (more than 10 cm in width) and not positioned inside a pressure vessel. Furthermore, the tomographic reconstruction algorithm assumed straight ray paths (Dines & Lytle 1979). In addition to the limitation imposed by the latter hypothesis when strong (more than 10 per cent) velocity contrasts develop, the method used by Falls et al. (1992) and Jansen et al. (1993) is not easily extendable to 3-D tomography, which would require a much larger number of sensors to provide a reasonable ray coverage of the samples. Due to space and sample size limitations, typical triaxial deformation apparatus are equipped with only 10-20 piezoelectric transducers, thus making purely active tomography impossible in practice (see Brantut 2010, section 8.9).
Here, I present a method that combines both active ultrasonics and passive AE monitoring data to improve the ray coverage in laboratory samples, and achieve jointly 3-D, anisotropic tomographic reconstruction and improved AE event location. During brittle failure or cataclastic flow experiments, the number of AE events is several orders of magnitude larger than the number of active transducers; furthermore, the AE locations are by definition more widely distributed (even though they often become clustered when strain localizes) than the fixed positions of the sensors. Therefore, the key idea is to use the AE events as passive sources in order to increase dramatically the ray coverage of the sample. Unfortunately, this major advantage comes with a significant cost: the determination of the AE source locations becomes coupled to the determination of the velocity structure. In other words, one cannot simply assume that the AE locations are known exactly, so that they become part of the inverse problem. Translated into the language of inverse problems, the observed data d are the arrival times measured from active and passive recordings, and the model parameters m are the wave velocities and anisotropy parameters at discretized points in space within the material, together with the source locations of all AE events. The physical model employed to make predictions of d based on known parameters m, is an eikonal solver g.
The principle of this approach is not new. It is essentially similar to Local Earthquake Tomography, a method developed in seismology to determine local velocity structures and accurate earthquake locations in seismically active regions (Aki & Lee 1976;Thurber 1983). However, the use and implementation of the technique in a laboratory setting presents several important challenges. The most prominent one is the efficient computation of theoretical arrival times (the forward problem) in three dimensions, including anisotropy and with the capability of handling strong contrasts. This combination of features is rarely required in seismological studies, and it was therefore deemed necessary to implement a specific method that could handle it. The subsequent key step is to com-pute the sensitivity of the arrival times with respect to the inversion parameters: velocity and anisotropy in the material, as well as source locations. Here again, a specific method using posterior ray tracing and trilinear interpolation was developed to perform this step efficiently. One other specifics of the laboratory studies is that the array of piezoelectric transducers is used for both passive and active measurements (whereby they act as ultrasonic sources as well as receivers). It is therefore important to include both types of data, and to accurately include their different range of measurement errors (i.e. variances and covariances on d and m). The very large number of AE events typically observed during brittle failure or cataclastic deformation of rocks (at least several thousand events) is a great advantage in terms of data coverage, but imposes the need to use automatic P-wave arrival time picking techniques which are never perfect. The challenge here is to achieve accurate picking of a majority of waveforms, while selectively rejecting bad picks to avoid including outliers with large errors in the inversion scheme.
Taken together, the specifics associated with using laboratory data call for the development of a specialized computational tool to perform what could be called 'Local Acoustic Emission Tomography'. In this paper, I present a set of methods and their implementation into an integrated free software called FaATSO (acronym for Fast Marching Acoustic Emission Tomography using Standard Optimisation) that performs 3-D, anisotropic AE tomography. The code is based on two key algorithms: (1) the fast marching method (Sethian 1999), employed to solve the forward problem, and (2) the quasi-Newton method (e.g. Tarantola 2005), used to solve the inverse problem. The software package is written in C++ with Matlab wrappers, and released online (under a free license, GNU GPL version 3) at https://github.com/nbrantut/Faatso.git. The performance and potential applications of the method are demonstrated on experimental data obtained during compactant deformation of a porous sandstone.
The paper is organized as follows. In Section 2, I describe the fast marching method and its implementation for a simple case of elliptical anisotropy relevant to most rock deformation experiments. Section 3 presents the quasi-Newton method to solve the inverse problem, and a particular emphasis is placed on the subtleties associated with the determination of ray paths and sensitivity kernels. A complete application of the method for a sandstone compaction experiment is given in Section 4. Conclusions and perspectives are presented in Section 5.

F O RWA R D P RO B L E M : 3 -D FA S T M A RC H I N G M E T H O D W I T H E L L I P T I C A L A N I S O T RO P Y
A wide variety of numerical methods exist to compute arrival times of seismic (or acoustic) waves in complex media (e.g. Rawlinson & Sambridge 2003). Because of the strong heterogeneities expected in laboratory samples, a natural choice is to use methods based on full wave front propagation rather than ray tracing. One of the most computationally efficient algorithms to compute wave fronts in strongly heterogeneous media is the fast marching method, which has been developed to solve so-called 'viscosity' solutions of Hamilton-Jacobi equations; a class of partial differential equations that includes the eikonal equation for seismic traveltimes (see review in Sethian 1999). The key benefit of the fast marching method is its unconditional stability, regardless of the velocity heterogeneity encountered by the wave front. Furthermore, wave front tracking is also beneficial because a single computation yields arrival times at many possible receivers from a given source. In this section, I Acoustic emission tomography 2179 briefly recall the principle of the method and explain the key implementation details for transversely isotropic media with vertical axis of symmetry (VTI).

Principle
For an anisotropic medium, the eikonal equation reads: where u is the wave front arrival time, x = (x, y, z) denotes Cartesian coordinates and V is the wave speed, which depends on the position and the orientation of the wave front. For VTI materials, the phase velocity is entirely determined by the phase angle θ , which is the angle between the direction normal to the wave front and the vertical direction. A discretized version of eq. (2) using upwind first-order finite differences on a regular Cartesian grid is due to Rouy & Tourin (1992): where V ijk (θ ) is the wave velocity at node (i, j, k), and denote first-order finite-difference operators, where h is the (uniform) grid spacing. The fast marching method consists in solving eq. (3) sequentially, from small to large values of u, hence progressively advancing the wave front. The use of upwind finite differences is absolutely key: it allows the preservation of causality and ensures the unconditional stability of the method (if first-order finite differences are used). A discussion and proof of why this holds is given in Sethian (1999). The method was originally developed for isotropic eikonal equations, but can also be applied for anisotropic equations without modification in special cases (further details below). The algorithm is the following: (i) Fix the starting source point at u = 0, and tag this point known, all other points begin tagged as unknown.
(ii) Update the neighbours of the newly known point by solving eq. (3), and tag them trial and remove them from the unknown set.
(iii) Choose the point with minimum u among the trial points, tag it as known and remove it from the trial set.
(iv) Go back to (ii), until no points remain tagged as trial.
The key step is (iii), which ensures that points are updated with increasing values of u (the information does not flow backwards). The computational efficiency of the method relies on the storage of the trial points in a so-called 'min-heap' data structure (in FaATSO, a Fibonacci heap from the C++ Boost library is used), so that the trial point with minimum arrival time is always at the top of the heap, and the addition of trial points and reordering of the heap is performed more efficiently than with conventional sorting algorithms.

Anisotropy
In case of VTI materials, the velocity is only dependent on the angle θ between the vertical axis and the direction of the normal to the wave front. In this work, I use a simple elliptical anisotropy: parametrized with the horizontal wave velocity, V 0 , and the relative excess in vertical wave velocity E. The parameter E is related to Thomsen's parameter : assuming a model for weak anisotropy and δ = (Thomsen 1986). In this model, V(θ ) is the phase velocity and θ is the phase angle. The cosine of θ is related to the partial derivatives of u in the following way: Although the fast marching method requires some modifications for general anisotropy (e.g. Sethian & Vladimirsky 2003), it can be shown that the conventional method as outlined above is valid when the direction of anisotropy is aligned with the grid and the anisotropy is not too strong. This is the case here provided that −1/2 < E < 1 (see Appendix A), which is largely in the remit of the approximation of weak anisotropy.
In practice, the assumption of VTI geometry is particularly relevant to initially isotropic rock samples subjected to triaxial compression in the laboratory. Indeed, loading under triaxial conditions tend to generate stress-and crack-induced anisotropies aligned with the compression axis (Paterson & Wong 2005, see also Section 4).

Second-order finite differences
The first-order finite-difference approximation of the eikonal equation is quite crude, and tends to generate large errors along diagonal points of the grid (e.g. Rickett & Fomel 1999). This can be mitigated by using a second-order (upwind) approximation, as described in Sethian (1999). The discretized eikonal equation becomes: where s i, j, k are switches defined by and One important fact to keep in mind when using a second-order approximation is that the scheme may not be applicable when sharp changes in wave speeds occur. In that case, the solution of u i, j, k might not exist (or would be a complex number), and we have to revert to the first-order scheme [eq. (3)].

Ray tracing
Once the arrival times have been computed at all points, it is relatively easy to compute ray paths a posteriori. In the isotropic case, ray tracing is straightforward because the ray angle at every point 2180 N. Brantut is perpendicular to the wave front (i.e. following the gradient of the arrival time field). In the anisotropic case, the situation is somewhat different because the rays are no longer perpendicular to the wave front: the phase angle is not the same as the group angle. So, at every point, one needs to compute the group angle (from the phase angle, given by the gradient of the arrival time field) in order to obtain the ray orientation. The relationship between group velocity V group and phase velocity V is (Thomsen 1986): where φ is the group angle and θ is the phase angle. The relationship between the group and the phase angles is: Here, the phase velocity is given by Therefore, the group angle is then given by and the group velocity is (as a function of the phase angle θ) where Fig. 1 provides an illustration of the differences between group and phase angles for an anisotropy parameter of E = 0.25. In practice, the rays are traced from the receivers to the source by successively generating piecewise linear segments oriented according to the local group angle as computed from eq. (14). The segment length is chosen to be equal to the grid spacing.

Accuracy
Strictly speaking, the finite-difference scheme is not globally second-order accurate, because the conditions for the switches are not always met. This is usually quite benign, except for the points immediately surrounding the source node: this is where the error is the largest, especially along the diagonals, because there is no way to use a second-order scheme around the source. Therefore, if used as is, the large first-order error generated around the source node will be 'dragged' along throughout the grid, even if secondorder approximations are made after that. This is illustrated in Figs 2 (a) and (b), which show a very large error of the fast marching method results (up to 20 per cent) along the diagonal points near the source region. A way around this issue can be found in Rickett & Fomel (1999), and consists of computing arrival times immediately around the source node analytically, assuming a constant velocity (and anisotropy) there. This 'source box' is made large enough to include all the eight nearest points within a cube surrounding the source, and especially includes all the first diagonal points. The improvement resulting from this technique is illustrated in Fig. 2(c), which shows a significant reduction of the error near the source region and throughout the medium (less than 5 per cent) when the source box is used compared to when it is not.
The overall accuracy of the numerical scheme is explored in Fig. 3, where a complex heterogeneous and anisotropic velocity Figure 2. Accuracy of the second-order finite-difference scheme used in the fast marching method for a homogeneous, anisotropic case. (a) Contours of computed arrival times using the analytical solution (black) and the fast marching method with (red) and without (blue) a source box. The horizontal wave velocity is V 0 = 1 mm µs −1 and the anisotropy parameter is E = 0.25. The medium is discretized in an array of 21 × 21 nodes. Relative errors are shown for the case without (b) and with (c) the source box. structure is used. Due to the strong heterogeneities, with V 0 ranging from 0.5 to 5.5 mm µs −1 and E ranging from −0.32 to +0.40, the wave front deviates from circularity and becomes cuspate in the low-velocity zone. The relative error in arrival time at receivers positioned on the boundary of the model (displayed as black dots) is computed as a function of grid spacing from 0.05 to 1 mm (i.e. from 401 to 21 gridpoints in each direction). In the absence of a true analytical solution for this velocity structure, a model with grid spacing of 0.02 mm (1001 points) is used for reference. The average error ranges from around 2 per cent in the coarsest grid, down to 0.03 per cent in the finest grid. The error decreases linearly with decreasing step size, indicating that the numerical method is globally of first order; this is not surprising considering that switches (eq. 8) are not always activated.
The computation times are remarkably short, from around 0.001 s in the 21 × 21 grid to 1.75 s in the 1001 × 1001 grid on a laptop using a single core of a 2.8 GHz Intel i7 processor. The speed of computation does not depend on the model complexity, but computations are moderately faster in the absence of anisotropy. As explained below, a grid made of around 1 million nodes is a con-venient choice in practice to achieve computations in a reasonable timescale.

Applicability and limitations
The computation of arrival times is performed by solving the eikonal equation, which is a high-frequency approximation of the wave equation. Therefore, there are several limitations to bear in mind regarding the applicability and physical resolution expected when using this technique to simulate laboratory observations.
In laboratory studies, AEs and ultrasonic data are often recorded using uncalibrated, resonant piezoelectric transducers the characteristics of which are given by the piezoelectric properties and dimensions of the transducer and casing. The resonant frequency of typical transducers is of the order of 1 MHz, which is characteristic of lead-zirconate-titanate ceramic discs of 1-2 mm in thickness and 4-5 mm in diameter. In a rock where the P-wave velocity is 5 mm µs −1 , the wavelength of the signals is of the order of 5 mm. The eikonal approximation is only valid if spatial variations in velocity occur over length scales significantly larger than the wavelength, which imposes a physical lower bound on the spatial resolution of the velocity model. Any wave velocity variations occurring over length scales smaller than the signal wavelength, that is, around 5 mm, will be averaged out, and may produce scattering.
In practice, AE and ultrasonic monitoring is typically performed using laboratory samples of 4-10 cm in diameter and 8-20 cm in length, so that the high-frequency approximation is often justified.
The extension of the method to finite frequencies would require the computation of the full wavefield to solve the forward problem, which is dramatically more computationally expensive than the computation of high-frequency wave fronts. Recent advances towards understanding waveform records from active measurements (e.g. Fukushima et al. 2009;Yoshimitsu et al. 2016) show that the radiation patterns from piezoelectric transducers (see Tang et al. 1994) and the cylindrical geometry introduce significant complexity in the transmitted signals and enhance the measurement errors. Such measurement errors, as well as the imperfection of the forward problem itself, can be accounted for in the formulation of the inverse problem (see below).

Assumptions and prior information
The fast marching method is an efficient tool to compute arrival times throughout a material based on a known velocity structure and source position. The data set used for tomography consists of arrival times from both active and passive sources. In typical laboratory experiments, the number of AE events (hereafter denoted n evt ) is much larger than the number of transducers (denoted n chan ). Taking advantage of the reciprocity of the wave equation, it is sufficient to compute arrival times using the transducers as sources and collect the arrival times at the positions of the events. The same set of transducers is often used as sources in the active measurements, so that there are up to n chan × (n chan − 1 + n evt ) arrival times. In practice, not all arrivals can be accurately picked on all channels, so that the actual quantity of data is somewhat reduced.
As stated previously, the model parameters are the values of V 0 and E throughout the material, and the source locations. The material is discretized into nodes along an orthonormal grid, with the z−axis being the vertical direction. The choice of the number of nodes is dictated by several constraints. From a physical point of view, the size of heterogeneities sensed by elastic waves is limited by their wavelength. In laboratory rock deformation experiments, the signal frequency is of the order of 0.5-1 MHz, so that the wavelength is of the order of several millimetres. Instead of using a relatively coarse discretization with independent nodes, here a much finer grid is used but with a correlation between neighbouring nodes over a fixed length scale. This length scale is at least commensurate to the wavelength of the signal used to probe the system, but can be adapted to the known type of heterogeneity of the material (Tarantola 2005). Furthermore, the use of a fine grid would also enhance the accuracy of the forward model computations. However, it is also desirable to limit the quantity of unknowns (especially since they are correlated) in order to ease the numerical computation and linear algebra involved in the determination of the solution (see the details below). Here, these contradictory requirements are met by using a relatively coarse grid to define the model parameters, and by refining it with trilinear interpolation when computing the forward problem. This approach allows us to minimize the number of parameters while retaining good accuracy in the computation of arrival times.
Typical cylindrical laboratory samples of 40 mm in diameter and 100 mm in length are discretized with an initial grid spacing of 5 mm and a refinement factor of 8, thus yielding around 2800 nodes (including additional grid steps outside the region of interest) in the coarse grid and more than 1.4 million nodes in the refined grid.
Due to the large number of unknowns (several thousands) and the computational cost of forward modeling, the inverse problem is not easily solved by grid search or Monte Carlo sampling methods, and we have to resort to least-squares optimization techniques. In this context, all observational errors and a priori information on model parameters are assumed to be described by Gaussian statistics.
The observed data (arrival times) vector is denoted d obs . The covariance matrix on the data, C D , is assumed to be diagonal (i.e. the observations are independent from each other). The diagonal elements of C D correspond to the variances on the observed arrival times. These variances are effectively the estimated errors on the observations. The passive and active measurements do not have the same errors, and two distinct variances are introduced: one for active survey data σ 2 shot , and one for the passive AE arrival time data σ 2 evt . The latter is typically larger than the former, since the active survey arrival times can be picked very accurately using cross-correlation methods on very clear waveforms, whereas AE arrival times are picked automatically and just manually checked. In general, it would be preferable to associate each arrival time with a specific error, for instance using the information provided by the autopicking algorithm. Such a method would require further developments that are beyond the scope of this work, so here I have used a single value of σ evt for all measurements.
In the framework of probabilistic inverse problems, the intrinsic error on the theory (here the ray approximation) can also be accounted for (Tarantola 2005, section 1.3). In practice, this is done by adding an additional variance on the data, that is, increasing the observational errors to include the possibility that systematic errors can be made due to the imperfection of the theory even if the observations are perfect.
The vector of a priori model parameters is denoted m prior . The prior information is also assumed to be Gaussian. Because the wave velocity is a strictly positive parameter [or a so-called Jeffreys parameter, see Tarantola (2005)], the logarithm of V 0 is used as the corresponding Cartesian parameter. All other parameters (anisotropy and source positions) are Cartesian. The covariance matrix in the model space C M is determined by: the variance σ 2 v and covariances on ln (V 0 ), the variance σ 2 E and covariances on E and the variances σ 2 x , σ 2 y , σ 2 z and σ 2 t0 on source locations (I assume no covariance between source location parameters). The covariances on ln (V 0 ) and E between two nodes of indices i and j follow a 3-D exponential covariance: where (x i , y i , z i ) are the coordinates of node i, σ 2 is the variance on either ln (V 0 ) or E and λ is a correlation length. Essentially, the length λ provides a smoothness of the velocity structure, preventing (statistically) heterogeneities over length scales smaller than λ.
In summary, the covariance matrix C M is block diagonal, the first block being the c ij for the logarithm of the velocity, the second block is the c ij for the anisotropy parameter and the remaining block is a simple diagonal containing the sequence of variances on source locations.
The inverse problem is non-linear: the function g(m) in eq. (1) effectively contains integrals of slownesses over ray paths between sources and receiver. The tomographic problem is only weakly nonlinear, while the source location problem is strongly non-linear. Because a Gaussian model is assumed for all parameters, it is essential that the a priori source locations are not too far from the true ones. In practice, initial source locations are determined using straight ray paths and an average P-wave velocity model (possibly anisotropic if needed).

Quasi-Newton algorithm
The method chosen to determine the mean posterior model parameters (i.e. the 'best' model in the sense of the least-squares norm) is the quasi-Newton method. It is an iterative quadratic minimization method, described in detail in Tarantola (2005, p. 78). The model parameters m are updated at each iteration n as where d n = g(m n ) and μ n ≈ 1 is the step size. The matrix G n contains the derivatives of g: In terms of implementation, the linear algebra operations in eq. (18) are performed using the C++ Armadillo package (Sanderson & Curtin 2016). Only the inverse of C M and C D are found explicitly, and a linear system is solved at every step instead of computing the inverse of G t n C −1 D G n + C −1 M . The key difficulty in the inversion procedure is computing G n .

Computing the matrix of derivatives
The (vector) function g maps the velocity structure and source locations to the arrival times at each receiver. If we denote t ij , the arrival time at receiver j from source i, the mapping d = g(m) can be written as where t 0 is the origin time of the source, R i j is the ray path between source and receiver, s is the curvilinear coordinate along the ray  path and V group is the group velocity along the ray. The medium we are dealing with is VTI, so that the group velocity depends on the local group angle φ between the ray and the vertical axis. The group velocity is expressed as where θ is the phase angle, easily computed and stored during the ray tracing process, and F is given by eq. (16). Denoting v = ln (V 0 ), Eq. (20) is rewritten as and its discretized form is then where the summation is performed over all the indices k = 0, . . . , N of the segments of length s k forming the ray R i j . The derivatives are therefore: and correspond to the orientation of the tangent to the ray at the source point, (x r0 , y r0 , z r0 ) being the coordinates of the first point along the ray originating from the source. As stated previously, the model parameters for the velocity structure are defined on a relatively coarse grid (Fig. 4a), which is refined by trilinear interpolation during the computation of the forward model (Fig. 4b). Rays are traced between sources and receivers in the refined grid, where ray segment lengths are uniform and equal to the (refined) grid spacing. The wave velocity, anisotropy, phase angle and ray point coordinates are stored, and the derivatives are computed with respect to the interpolated parameters. In order to access the derivatives in the initial coarse grid, which correspond to the elements of G n , the derivatives obtained on the refined grid are combined using the linear interpolation coefficients used in the grid interpolation procedure (Fig. 4c). For any model parameter on the coarse grid m j (either the horizontal wave velocity or the anisotropy parameter), we have where p k are the model parameters on the refined grid obtained from trilinear interpolation: where a kj are the interpolation coefficients, so that and we obtain

Synthetic examples
In order to assess the efficiency and accuracy of the inversion method, it is first used on synthetic data sets for which the true model parameters are known. The synthetic data are generated from a 2-D heterogeneous velocity field with an average of V 0 = 3 mm µs −1 and up to 15 per cent spatial variation (from around 2.8 to 3.2 mm µs −1 ). The full grid size is 5 × 5 mm 2 , and the velocity is  Table 1. set to 0 outside a circular domain of 5 mm in diameter. The nominal grid spacing is 0.2 mm, and a refinement factor of 10 (grid spacing of 0.02 mm) is used in the forward computations. A set of 40 AE events are randomly positioned in the material, and eight sensors are located on the circular edge of the material. The exact arrival times are computed with the fast marching method, and a random error (±5 per cent) is added to simulate measurement uncertainty. The prior velocity structure is assumed constant, and the prior AE event locations are given by the true ones plus a random error of ±0.2 mm. A first example is computed in a purely isotropic case, as shown in Fig. 5. The parameters used in the inversion are summarized in Table 1, and the homogeneous prior model is shown in Fig. 5(a). The quasi-Newton algorithm converges after five iterations, and least-squares residual drops by a factor of around 23 between the initial and final iterations. The result of the inversion is the 'best' solution in the sense of the least-squares norm, shown in Fig. 5(c), and corresponds to the average of all likely models once the information contained in the observations is used (Tarantola 2005). The comparison of the true model (panel b) and inverted model (panel c) reveals that the inversion procedure allows us to retrieve the main features of the true model, even though the prior model is uniform. The roughness of the inversion result is controlled by the combination of the variance σ v and the correlation length λ. Increasing λ or decreasing σ v reduces both the amplitude of the changes in V 0 and the length scale over which the changes are observed. In all cases, the AE event locations are retrieved almost exactly after a few iterations, despite erroneous initial locations. The average relative error between the true and inverted locations is around 6 per cent. The improvement in the AE event locations is quite modest, and the inversion method stops converging if the a priori locations become too different from the true ones. This is a clear limitation of the method and limits its applicability regarding AE locations alone. However, the key benefit provided by the inversion is that we can retrieve the main features of the velocity structure.
Another example is shown in Fig. 6, where heterogeneities in both V 0 and E are implemented. The same geometrical pattern is used for velocity and anisotropy. The parameters used in the inversion are summarized in Table 1, and the homogeneous prior models are shown in Figs 6(a) and (d). Similarly to the isotropic case, the quasi-Newton algorithm converges after five iterations, and least-squares residuals drop by a factor of around 93 between the initial and final iterations. The average models resulting from the inversion procedure are shown in Figs 6(c) and (f). Again, the main features of the true model are recovered by the inversion. Because now both V 0 and E are inverted, it is important to determine whether the inversion results are correlated; in other words, it might be that jointly changing V 0 and E allows us to fit the observation equally well, in which case there would be a trade-off between these parameters. In order to estimate the posterior correlations between model parameters, we compute the posterior covariance matrix as (Tarantola 2005) and the correlation matrix as The elements of ρ corresponding to the correlations between the velocity V 0 and the anisotropy E at the same node are shown graphically in Fig. 7. All the correlations are negative, and their magnitude is relatively large (down to −0.5) where the ray coverage is poor, typically between pairs of sensors where no AE events occur. Overall, the correlations are quite small (−0.1 to −0.2), where the number of AE events is significant, such as at the centre of the model. The synthetic tests shown here are not intended to cover all possible situations and there are known caveats associated with this exercise (Tarantola 2005;Igel 2017). Nevertheless, the overall very good match between the inversion results and the known true models gives us confidence that the method is able to retrieve the main  features of the velocity structure based on AE and active measurements. The output (mean) models are spatially smooth, a feature that is entirely controlled by the variances σ v and σ E and the correlation length λ. The true models used here are also smooth, and the good performance of the inversion method strongly depends on the a priori knowledge of that smoothness, so that the variances and correlation length can be tuned to the expected variations. In practice, when working with real data, this tuning is left to the user, and any useful prior knowledge (including: signal wavelength, sample geometry, macro-and microstructural observations and initial AE event locations) should be used to constrain the covariance matrix.

A P P L I C AT I O N T O C O M PA C T I O N O F P O RO U S S A N D S T O N E
In order to test the tomographic method on real data and demonstrate its potential for laboratory rock deformation experiments, I perform a sequential inversion of AE and ultrasonic data obtained during a triaxial deformation experiment conducted on a porous sandstone. Sandstone is probably the ideal rock type for such an analysis because of the very large number of AE events and strong variations in wave velocity typically observed during its deformation (e.g. Fortin et al. 2006). Furthermore, since ultrasonic measurements can be performed repeatedly during a single experiment, if enough AE events are recorded during the time window around each active measurement, the inversion is expected to provide a complete sequence of wave velocity evolution throughout the deformation process.
In porous sandstones, the mode of deformation transitions from brittle to ductile with increasing confining pressure. Ductile deformation is initially compactant and proceeds by cataclastic flow, which is driven by microscale fracturing. While brittle deformation is always associated with localized faulting, the occurrence of strain localization during ductile deformation depends on both the initial microstructure of the rock and the stress state at which the rock is deformed (Paterson & Wong 2005;Wong & Baud 2012). In particular, a regime of compaction localization is observed at elevated confining pressures in some high porosity sandstones (e.g. Klein et al. 2001;Baud et al. 2004Baud et al. , 2006. The formation of compaction bands in sandstone was analysed in detail by Fortin et al. (2006) using AE monitoring and ultrasonic measurements, which revealed that the growth of the bands coincides with clusters of AE events and is accompanied by a progressive drop in the average wave velocity in the material. Here, I performed a triaxial deformation experiment on a sandstone similar to that used by Fortin et al. (2006), in the compaction banding regime. As explained in detail below, the use of AE tomography in this experiment allows us to resolve the local changes in wave velocity associated with the propagation of compaction bands, which was previously beyond our capabilities.

Experimental method
The rock chosen for this study is a 22 per cent porosity sandstone from Bleurville (Vosges, France). This material is composed of a majority of quartz and K-feldspar, and comes from the same area as the Bleurswiller sandstone used by Fortin et al. (2006). A cylinder of 40 mm in diameter was cored perpendicular to bedding, and was then cut and precision ground to 100 mm in length. The sample was inserted in a viton jacket and positioned between two steel end caps within the triaxial rock physics apparatus of the Rock and Ice Physics Laboratory at University College London (see description in Eccles et al. 2005).
The sample was saturated with distilled water. The confining pressure and pore pressure were set to 90 and 10 MPa, respectively, and maintained constant throughout the test by the use of servo-controlled intensifiers. The Terzaghi effective pressure was therefore equal to 80 MPa, which is well within the compactant regime for Bleurville sandstone. The pore volume change in the sample was monitored by measuring the change of fluid volume in the pore pressure intensifier with an LVDT (Linear Variable Differential Transformer) attached to the moving actuator. The axial load was measured with an external load cell and the sample shortening was measured with an external LVDT, corrected for the elastic compliance of the loading column and sample assembly. The deformation test was run at a constant strain rate of 10 −5 s −1 .
The viton jacket positioned around the sample was equipped with an array of 16 piezoelectric transducers arranged in the geometry shown in Fig. 8. The piezoelectric elements of the transducers are Figure 9. Differential stress (black), porosity reduction (grey) and AE rate (green) as a function of axial strain during triaxial deformation of Bleurville sandstone at 80 MPa effective confining pressure. Black diamonds indicate the points where active ultrasonic measurements were performed. discs of lead-zirconate-titanate ceramic (material reference PIc255 purchased from Physik Instrumente GmbH) of 2 mm in thickness and 5 mm in diameter. These discs are mounted on top of 12 mm long aluminium inserts which act as waveguides; the additional traveltime within the waveguides was calibrated using a plain aluminium rod of known properties as the sample.
The signals from the transducers were amplified by 40 dB and recorded at 50 MHz with a digital oscilloscope (12-bit dynamic range). Passive AE monitoring was performed by recording and storing 4096 point waveforms on all channels each time a signal with an amplitude higher than 120 mV was detected on at least two channels, up to a maximum of around four events per second. The background noise level was around 20 mV. In addition, repeatedly during the experiment, active ultrasonic monitoring of the sample was performed by sending a high-frequency (1 MHz) high voltage (250 V) pulse to each transducer (then used as a source), and using the other sensors in the array to record the waveforms. Each shot was performed six times and the resulting waveforms were stacked to enhance the signal-to-noise ratio.
The waveform and amplitude spectrums of a typical AE event is shown in Fig. B1, where we observe that most of the signal is at a frequency of around 0.4 MHz. The P-wave velocity of the sample is of the order of 4 mm µs −1 , so that the wavelength of the signal is of around 1 cm. This value is four times lower than the sample diameter, and 10 times lower than the sample length.

Data analysis
The resulting mechanical data are presented in Fig. 9. The differential stress on the sample (Q) initially increases linearly with increasing axial strain, while the porosity also decreases. The onset of inelastic compaction is at around 45 MPa differential stress (i.e. 0.36 per cent axial strain). Beyond that point, the porosity decreases non-linearly with increasing stress. At around 1 per cent axial strain, the differential stress deviates significantly from linearity and stabilizes at around 105 MPa. The stress remains approximately constant between 1.2 and 2.2 per cent axial strain, and linear strain hardening is observed beyond that point, with the stress reaching around 120 MPa at 5.5 per cent axial strain. This mechanical behaviour is typical of the ductile, compactant regime in porous sandstone (Wong & Baud 2012). A total of 10 163 AE events were recorded during deformation. The AE rate as a function of axial strain is shown in green in Fig. 9. Based on their high amplitude and the clarity of the first arrivals, a subset of 97 740 waveforms from 6458 AE events was selected to be used in the tomographic inversion. The arrival times were initially picked automatically using a method based on an autoregressive time-series model, in which the number of autoregression coefficients was determined by the Akaike Information Criterion (Sleeman & van Eck 1999). All picks were then individually checked and corrected manually.
In addition to the AE events, active ultrasonic surveys were performed at the times marked by black diamonds in Fig. 9. The waveforms associated with an initial survey prior to deformation were picked manually, and the changes in arrival time for all subsequent surveys were computed automatically using waveform crosscorrelation with subsampling at 100 MHz (Brantut et al. 2011;Brantut 2015). This procedure yielded a set of 147 accurate relative picks for each survey (with redundant information considering that each transducer is used as both a source and receiver sequentially).
The time-series of AE and survey data is split into 16 sequences centred around the time of each survey (labeled 'a' to 'p' in Fig. 9), leading to 16 individual data sets each containing 1 survey and from 50 (step 'a') to 625 (step 'c') AE events. The first two surveys (before 'a') were not used due to the lack of high-quality AE events. A tomographic inversion is performed for each of these data sets. The parameters used in the inversion method are summarized in Table 2. The inversion of the first data set ('a') is performed by assuming an a priori velocity structure that is homogeneous, anisotropic with V 0 = 3.9 mm µs −1 and E = 0.01. These values correspond to the average velocity model as measured by the active ultrasonic survey. The inversion of data sets 'b' to 'p' was then performed sequentially, by using the output velocity structure of the previous set as the a priori model for the current set. The a priori covariance matrix was left unchanged. In each set, the prior AE locations were set to those found with a grid search approach (with least-absolute value minimization) in the homogeneous, anisotropic velocity model as determined with the ultrasonic surveys.

Inversion results
The results of the sequential inversion procedure are summarized in Figs 10 and 11. Fig. 10 shows the spatial evolution of the relative change in horizontal velocity, computed as (V 0 − V (a) 0 )/V (a) 0 , in a vertical slice cutting through the centre of the sample along the x-direction. AE hypocentre locations for each data subset are also shown on each subplot. During step (b), AEs appear quite scattered throughout the upper part of the sample, and the velocity structure is relatively homogeneous, except from a very slight anisotropy also localized at the top of the sample. By contrast, during step (c), the AE locations form a clear subhorizontal band located near the top of the sample. This clustering of AE hypocentres is accompanied with a drop of the wave velocity and an increase in anisotropy around the AE cluster. During each subsequent step from (d) to (n), the AE hypocentres form a subhorizontal front that gradually propagates from the top to the bottom of the sample. Behind the AE front, the wave velocity drops quite abruptly by 5-15 per cent and the anisotropy increases by 5-15 per cent as well, indicating that the vertical P-wave velocity remains almost unchanged, while the horizontal wave speed is strongly affected. During steps (o) and (p), the AEs appear less clustered and the wave velocity structure becomes more homogeneous, albeit with some scattered patches of slightly higher velocity.
As explained from the synthetic tests, the overall smoothness of the velocity structure is controlled by the combination of the variances σ v and σ E and the correlation length λ. The choice of a rather large λ (50 mm) was made to overconstrain the model and ensure relatively small changes between each step. Despite imposing this level of smoothness, relatively sharp transitions are observed, for instance during steps (g) and (h). This is because the inversion is performed sequentially, using the successive outputs as a priori models for the following steps.
Obviously, changing the a priori variances or correlation length would modify quantitatively the outcome of the inversion. An exploration of all possible outputs is beyond the scope of this work, but a number of tests were performed with a range of parameters (σ v from 0.12 to 0.3, σ E from 0.02 to 0.1 and λ from 20 to 200 mm), revealing that the pattern observed here, that is, the propagation of a front of AEs behind which the velocity drops, is retrieved systematically.
The reality of the patches of higher velocity imaged in the centre of the sample in steps (o) and (p) is difficult to confirm. The ray coverage is very poor in these areas, since all the AEs are located at the bottom of the sample and the sensors are on the sides. Hence, the active surveys do not provide enough coverage on their own to robustly determine the structure in the central part of the sample.

Significance
The combined AE locations and tomography produce an unprecedented insight into the mechanics of compaction in porous sandstone. The AE hypocentres delineate a succession of compaction bands starting at one end of the sample and systematically progressing towards the other end. Behind the compaction front, a small number of diffuse AEs are recorded. However, the tomography reveals that the material has been intensely damaged and cracked, as evidenced by the inferred large drop in horizontal wave velocity (up to 15 per cent). Here, the AEs occurring at the front have been used to image the whole sample, even in regions where little or no AE activity was recorded.
More diffuse AE activity is recorded in the last steps of deformation (o) and (p), and the tomography seems to produce a more homogeneous (anisotropic) velocity structure. The structure inferred from the last deformation step compares well with the macroscopic 2188 N. Brantut The reference velocity field is taken at step (a). All slices correspond to north-south cross-sections ((x, z) plane) running through the centre of the sample. Each slice is 40 mm in width and 100 mm in height, which matches the initial sample dimensions. The black dots correspond to the AE event hypocentres located within |y| ≤ 2.5 mm of the central section of the sample.
picture of the deformed sample recovered after unloading and decompression, as well as with the total number of AEs recorded in the sample, shown in Fig. 12.
These results imply that the accumulated AE density is correlated with the horizontal P-wave velocity drop. This is consistent with the idea that each AE corresponds to the propagation of one or several subvertical microcracks, which in turn are known to dramatically affect the wave velocity of the material (e.g. Sayers & Kachanov 1995). The correlation between the local cumulative AE number and the change in wave velocity is shown graphically in Fig. 13. There is a clear proportionality between the maximum AE number in a given location and the corresponding velocity change at that location. However, for a given velocity change, the cumulative AE number can be significantly less than the potential maximum, as evidenced by the spreading of the points along the y-axis in Fig. 13. The lack of direct proportionality between wave velocity change and AE Figure 11. Evolution of the anisotropy parameter E during deformation, as inferred from the inversion of data sets (a)-(p). All slices correspond to north-south cross-sections ((x, z) plane) running through the centre of the sample. Each slice is 40 mm in width and 100 mm in height, which matches the initial sample dimensions. The black dots correspond to the AE event hypocentres located within |y| ≤ 2.5 mm of the central section of the sample. count-the correlation being only with the maximum AE countis maybe an artefact arising from the AE detection and selection thresholds, but clearly highlights the potential of the tomographic method to image areas where no AEs are recorded but are still heavily microcracked.
The progression of the compaction front imaged by AE tomography is a key element to understand the mechanical behaviour of the rock. Throughout the propagation of the front, from steps (b) to (h), the applied differential stress remains constant (Fig. 9). The front seems to stagnate from steps (h) to (l), which corresponds to a period of slight strain hardening. Beyond step (l), the front seems to propagate further while losing its coherence, which corresponds to a decrease in strain hardening. During this stage, the whole sample becomes cracked and the wave velocity becomes again more homogeneous. During steps (h)-(l), the rock is effectively a bimaterial, with the top half heavily cracked and mostly incohesive, and  the bottom half more or less intact. The period of strain hardening indicates that the deformation of the top half requires progressively more stress to proceed [strain hardening is a common observation in the ductile cataclastic flow regime, see Wong et al. (1997)]. The subsequent progression of the compaction front in the bottom half of the sample, together with the slight decrease in strain hardening, shows that it becomes more energetically favourable to start compacting the intact part of the rock rather than continue deforming part the damaged part of the sample.
Overall, the tomographic imaging reveals a close relationship between AE locations and density, changes in P-wave velocity and compaction localization in the sample. In a previous study by Fortin et al. (2006) on a similar rock, AE locations were used to image the growth of compaction bands and the relationships between compaction bands and pre-existing heterogeneities in the sample. In that study, compaction bands initiated at various locations within the sample and did not start specifically at one end of the sample.
By contrast, here we find that compaction bands successively form next to each other from the top to the bottom of the sample. This key difference in deformation mode might be associated with the lack of pre-existing heterogeneities in our sample compared to that used by Fortin et al. (2006). This explanation is consistent with the initiation of compaction near the end of the sample, where the stress state is slightly more compressive due to frictional effects between the sample and the end-caps: this stress heterogeneity is the dominant one compared to other microstructural heterogeneities in the sample.

C O N C L U S I O N S
In this paper, I have presented a method that combines active ultrasonic monitoring and passive AE recording in laboratory samples to perform 3-D, anisotropic wave velocity tomography. The method is analogous to Local Earthquake Tomography, pioneered by Aki & Lee (1976) and Thurber (1983). The data used are arrival times from both active sources with known positions and passive sources with unknown positions. The model parameters are the wave velocities at regular gridpoints within the sample (a horizontal P-wave speed V 0 and an anisotropy parameter E, in an elliptic VTI model) and AE event locations. The forward problem is solved efficiently by the fast marching method (e.g. Sethian 1999). The inverse problem is solved with the quasi-Newton method (Tarantola 2005).
Tests performed using synthetic data with noise indicate a good performance of the method, with rapid convergence and good agreement between the inversion results and the true input model. The results are somewhat less accurate when performing inversion with anisotropy, and there are small negative correlations between the inverted wave velocity and anisotropy parameters. These problems are inherent to all tomographic methods and can only be circumvented by using better constraints on prior models and covariances.
The application of the method to experimental data obtained during compaction of a porous sandstone reveal unprecedented information on the relationships between strain localization, AE density and changes in wave velocity. The tomographic imaging shows very clearly the propagation of a compaction front, highlighted by clusters of AE events, and leaving a heavily damaged material behind (with drops in horizontal P-wave velocity of up to 15 per cent). These results are consistent with the post-mortem structure of the sample.
The largest source of error in the approach is, by far, the lack of accuracy of the P-wave arrival picks on passive AE records. Here, all the picks were manually checked, but this task would be near impossible on very large data sets (or, rather, not considered a good use of people's time). One way to improve this crucial step would be to use machine learning techniques to automate the picking and provide individual measures of picking errors for each arrival time. These individual errors could then be used to fill the data covariance matrix and obtain more robust inversion results.
The overall good performance of the method on real laboratory data demonstrates the potential of AE tomography as a routine in situ imaging technique. Indeed, the method requires only arrival time data from AE transducers, which are already routinely collected in several laboratories. It is therefore imaginable to revisit pre-existing data sets and produce new complementary results with AE tomography. To this end, I provide the full software suite FaATSO online with the required documentation and manuals (https://github.com/nbrantut/FaATSO.git), with the hope of encouraging the development and use of improved AE tomography techniques.

A C K N O W L E D G E M E N T S
The premises of this work were initiated a long time ago, and I am very grateful to Albert Tarantola for his teaching and his help in the very early stages of this endeavour. I also thank Romain Jolivet for key discussions on optimization methods, Ana Ferreira for her support and general advice on seismic tomography and Phil Meredith for his support and useful comments on an early version of this manuscript. This work was supported by the UK Natural Environment Research Council (grants NE/K009656/1 and NE/M016471/1). The software suite FaATSO can be accessed at https://github.com/nbrantut/FaATSO.git.