Imaging of Small‐Scale Heterogeneity and Absorption Using Adjoint Envelope Tomography: Results From Laboratory Experiments

To complement the information provided by deterministic seismic imaging at length scales above a certain resolution limit we present the first application of adjoint envelope tomography (AET) to experimental data. AET uses the full envelopes of seismic records including scattered coda waves to obtain information about the distribution of absorption and small‐scale heterogeneity which provide complementary information about the investigated medium. Being below the resolution limit this small‐scale structure cannot be resolved by conventional tomography but still affects wave propagation by attenuating ballistic waves and generating scattered waves. Using ultrasound data from embedded sensors in a meter‐sized concrete specimen we image the distribution of absorption and heterogeneity expressed by the intrinsic quality factor Q−1 and the fluctuation strength ɛ that characterizes the strength of the heterogeneity. The forward problem is solved by modeling the 2‐D multiple nonisotropic scattering in an acoustic medium with spatially variable heterogeneity and attenuation using the Monte‐Carlo method. Gradients for the model updates are obtained by convolution with the back‐propagated envelope misfit using the adjoint formalism in analogy to full waveform inversion. We use a late coda time window to invert for absorption and an earlier time window to infer the distribution of heterogeneity. The results successfully locate an area of salt concrete with increased scattering and concentric anomalies of intrinsic attenuation. The resolution test shows that the recovered anomalies constitute reasonable representations of internal structure of the specimen.


Introduction
Imaging methods that infer the internal structure of an object from measurements performed from the exterior are referred to as tomographic methods which are characterized by (a) the material property that is being imaged, (b) the physical observable that is measured and (c) the model for the physical interaction that connects (a) to (b). Examples for tomographic methods in geophysics include muon tomography (Lechmann et al., 2021) that uses the interaction of cosmic ray muons with matter to image the distribution of density or electrical resistivity tomography (Daily et al., 2004) that infers the subsurface resistivity structure from measurements of the electric field.
Seismic methods use the interaction of elastic waves with the subsurface to infer the distribution of elastic parameters. In the simplest case, the wave propagation velocity is reconstructed from observations of direct wave travel times (Aki, 1969). Depending on the interaction of elastic waves with the medium, an analysis of reflected waves is more successful-the approach of seismic reflection imaging (Gray et al., 2001) used to derive the subsurface impedance structure. A requirement for this approach to work is that the reflected arrivals are more or less Abstract To complement the information provided by deterministic seismic imaging at length scales above a certain resolution limit we present the first application of adjoint envelope tomography (AET) to experimental data. AET uses the full envelopes of seismic records including scattered coda waves to obtain information about the distribution of absorption and small-scale heterogeneity which provide complementary information about the investigated medium. Being below the resolution limit this small-scale structure cannot be resolved by conventional tomography but still affects wave propagation by attenuating ballistic waves and generating scattered waves. Using ultrasound data from embedded sensors in a meter-sized concrete specimen we image the distribution of absorption and heterogeneity expressed by the intrinsic quality factor Q −1 and the fluctuation strength ɛ that characterizes the strength of the heterogeneity. The forward problem is solved by modeling the 2-D multiple nonisotropic scattering in an acoustic medium with spatially variable heterogeneity and attenuation using the Monte-Carlo method. Gradients for the model updates are obtained by convolution with the back-propagated envelope misfit using the adjoint formalism in analogy to full waveform inversion. We use a late coda time window to invert for absorption and an earlier time window to infer the distribution of heterogeneity. The results successfully locate an area of salt concrete with increased scattering and concentric anomalies of intrinsic attenuation. The resolution test shows that the recovered anomalies constitute reasonable representations of internal structure of the specimen.

Plain Language Summary
No matter how small the structures are that a seismic imaging method is able to resolve, there is structure with smaller length scale. On the one hand this small-scale structure causes unwanted signals for conventional imaging approaches. But on the other hand it provides complementary information about the investigated medium. To turn this to our advantage we, for the first time, apply a new imaging method that uses the waves which are caused by the small-scale structure. Using data of an experiment in a concrete block we demonstrate that we can identify areas of anomalous small-scale structure. The results may help in the future to locate minute perturbations in the medium as they occur in the advent of volcanic eruptions or after earthquakes and to obtain new information about the geologic history of subsurface materials. The approach can be transferred to investigate man-made materials and structures, such as deteriorating concrete constructions.

ZHANG ET AL.
It is important to image the distribution of heterogeneity because it influences the wavefield by scattering which attenuates direct waves and generates coda waves (Aki, 1969;Aki & Chouet, 1975) which arrive at later times. In a more complex velocity structure as in the deep Earth, scattering can also cause precursory arrivals (Shearer, 2007). Besides its effect on the wavefield, the distribution of heterogeneity itself provides complementary information about the target medium which has been widely used in imaging of volcanoes (De Siena et al., 2013 and the deep Earth (Margerin & Nolet, 2003;Sens-Schönfelder et al., 2021). Since the value of tomographic imaging lies in the geological interpretation, the complementary information about heterogeneity can augment the conventional imaging of macroscopically averaged elastic properties.
Additional need for a tomographic method to image heterogeneity comes from the monitoring of elastic properties. Scattered coda waves have been shown to be highly sensitive to subtle changes of elastic properties which promoted many useful methods, such as coda wave interferometry (CWI, Poupinet et al., 1984;Snieder et al., 2002). Combined with seismic interferometry of ambient noise, this has been used for continuous monitoring of subtle changes in volcanoes (Sens-Schönfelder & Wegler, 2006), fault zones (Brenguier et al., 2008), environmentally stressed areas (Sens-Schönfelder & Eulenfeld, 2019), and even assessing groundwater storage (Illien et al., 2021;Mao et al., 2022). The spatial sensitivity of this coda wave-based monitoring, however, depends on the distribution of the heterogeneity that generates the coda waves (T. Zhang et al., 2021;van Dinther et al., 2021).
Meanwhile, the success of Snieder et al. (2002) inspired many laboratory experiments using ultrasonic waves in the concrete for monitoring the changes in materials corresponding to stress (Larose & Hall, 2009;Stähler et al., 2011), temperature (Larose et al., 2006;Niederleithinger & Wunderlich, 2013) or damage (Schurr et al., 2011;Wang et al., 2020). Planès and Larose (2013) reviewed the applications of ultrasonic CWI in concrete in detail. The development of CWI-related techniques benefited a lot from the parallel applications in acoustics, engineering and seismology.
A commonly used approach to describe wave scattering in heterogeneous media is diffusion theory that treats the propagation of wave energy analogous to the diffusion of for example, heat. In the diffusion theory, the diffusion constant and dissipation are used to describe the ultrasonic scattering and intrinsic attenuation (Anugonda et al., 2001;Becker et al., 2003). These two parameters can be estimated by comparison between the experimental data and theoretical predictions using the diffusion model. This allows people to describe the effect of uniformly distributed material damage (Deroo et al., 2010;Ramamoorthy et al., 2004). The diffusion model has also been used for calculating the sensitivity kernel of CWI to velocity changes or decorrelation which allows for imaging of the spatial distribution of the changes (Rossetto et al., 2011;Y. Zhang et al., 2016). Although the diffusion model has been successfully implemented to simulate the wave scattering and absorption, it is a simplification of the multiple-scattering process and hard to extend to more realistic cases, like the early coda, short source-receiver distances, anisotropic scattering or spatially variable heterogeneity. Wu (1985) first proposed the multiple scattering model and introduced the radiative transfer theory (RTT) to seismology. To numerically solve the radiative transfer equations, the Monte Carlo method was introduced (Gusev & Abubakirov, 1987;Hoshiba, 1991), which allows for the ability to simulate wave scattering in the spatially variable heterogeneity and intrinsic attenuating media (T. Zhang et al., 2021). Instead of the diffusion constant and dissipation used in the diffusion model, the spatial distribution of fluctuation strength ɛ in the random medium and the intrinsic quality factor Q −1 describe the spatial variability of scattering and absorption. The simulation of energy propagation with spatially variable properties using RTT allowed us to introduce the adjoint method initially developed in FWI (Fichtner et al., 2006(Fichtner et al., , 2010Tarantola, 1984;Tromp et al., 2005) for the imaging of scattering and absorption properties with scattered waves (T. Zhang & Sens-Schönfelder, 2022 In this article, we apply the new seismic tomography method that was introduced by T. Zhang and Sens-Schönfelder (2022) as adjoint envelope tomography (AET) to a real world data set. For this first application we selected a setting that (a) can be treated in two dimensions to reduce the computational demands in this first test, (b) allows for a treatment of scattering in the acoustic approximation and (c) provides independent information about the internal distribution of heterogeneity in the target region. We chose a metric-sized concrete slab with embedded ultrasound transducers and known internal structure as presented in Section 3. In Section 5, we discuss the significance of this engineering-scale experiment for applications to the Earth which is ensured by the common size-wavelength ratio of the structures imaged in the ultrasound experiment in concrete and solid Earth targets investigated with seismic waves.
Section 2 briefly introduces the methodology of AET. The experiment is described in Section 3 including the data processing and the investigation of the background parameters used as starting models in the iterative inversion. The inversions for absorption and scattering properties are conducted individually with later and early coda waves, respectively, in Section 4. Section 5 contains the interpretation of the inversion results, analysis of the resolution tests and a brief outlook for seismic applications.

Radiative Transfer Theory
Adjoint envelope tomography is based on RTT that describes the propagation of energy in scattering media. RTT was originally developed to investigate the propagation of light through atmosphere (Apresyan et al., 1996;Chandrasekhar, 1960) and was later introduced in seismology by Wu (1985). Assuming that the phase of interfering scattered waves is randomized, RTT uses the additivity of wave energy rather than wave amplitudes which allows to use a statistical description of the medium heterogeneity instead of a deterministic one that is required for the wave equation. For the goal of imaging the heterogeneity, it is necessary to model energy propagation in the presence of spatial variability in heterogeneity and intrinsic attenuation. The 2-D acoustic radiative transfer equation with spatially variable heterogeneity and absorption is written as (T. Zhang & Sens-Schönfelder, 2022): where ( , , ) is the specific energy density propagating in direction n which is part of the total energy density ( , ) at the position r with lapse time t. The left-hand side of Equation 1 is the derivative of the specific energy density in along a path element. The right-hand side describes the two processes that change the energy. The first term accounts for the loss of energy due to scattering from the given direction n into other directions (α 0 g 0 (ɛ 2 (r))) and due to intrinsic attenuation (ω/Q(r)). The second term describes the associated effect of scattering: the energy increase by scattering from all other directions into direction n. ( , ) can be compared to seismogram envelopes and constitutes the observable in AET. α 0 indicates the velocity of the acoustic wave and ω is the angular frequency. Q −1 (r) is the inverse intrinsic quality factor. g 0 (ɛ 2 (r)) is the total scattering coefficient that is defined as the angular integral of the scattering coefficient ( , ′ , ( ) ) as is given by Wegler et al. (2006): where k 0 is the wavenumber and θ is the scattering angle between the incident direction n and the direction of the scattered wave n′. Φ is the local power spectral density function (PSDF) of the spatial parameter fluctuations in the heterogeneous medium. The PSDF is the statistical characterization of the small-scale medium heterogeneity which influences the energy propagation through its effect on the scattering coefficient. In the present case an exponential type PSDF is assumed: with the wave vector k. As shown in Equation 3, the RTT uses two parameters to describe the heterogeneity of the propagation medium-the correlation length a and the strength of the parameter fluctuations ɛ. We assume 4 of 21 here that the correlation length a is constant throughout space and that the spatial variability of scattering and attenuation is fully described by the distribution of ɛ 2 (r) and Q −1 (r). This representation allows us to model the energy propagation in any spatially variable model = { 2 ( ), −1 ( ) } .

Adjoint Method and Iterative Inversion
The concept of AET is in full analogy to the adjoint method in waveform tomography (Fichtner et al., 2006;Tarantola, 1984;Tromp et al., 2005). The inversion of seismogram envelopes for heterogeneity and absorption models starts with the definition of the misfit function that quantifies how well the model predictions match the observed data. In the present case the misfit function measures the match between the observed envelope of ultrasonic waves D(r j , t; r i ) and the synthetic energy density E(r j , t; r i , m) simulated in the current model m. r i and r j represent the positions of the ith source and the jth receiver, respectively. Note that the envelope defined here is the squared velocity envelope. The least-squares misfit function is defined as: The integration time window [T 1 , T 2 ] can contain the ballistic waves or coda waves or it can comprise the full envelope including both ballistic and coda waves. But for simplification of derivation, we redefine the integration bound in Equation 4 as [0, ] where T = T 2 − T 1 . The goal of the inversion is to minimize χ(m) by adapting the model parameters m. This process requires knowledge of the Fréchet derivative that we denote δχ(m) and which describes the derivative of the misfit function χ with respect to changes in the model m. δχ(m) could be calculated explicitly using finite differences on each parameter in the model vector m separately, which, however, is very ineffective since the dimension of m usually is large. Instead we use the adjoint method that greatly simplifies the calculation of the Fréchet derivative. We first write the Fréchet derivative as an integral over all model parameters, that is, and integral over space V since we face an imaging problem: here ɛ K χ (r′) and Q K χ (r′) are the scattering and absorption misfit kernels, respectively, with respect to the changes in scattering and absorption properties δɛ 2 (r′) and δQ −1 (r′). In T. Zhang and Sens-Schönfelder (2022), we used the adjoint formalism to derive expressions for the misfit kernels: is the total scattering coefficient under the condition that ɛ = ɛ 0 , which involves normalization by 2 0 referring to Equations 2 and 3. f(n, n′) indicates the normalized differential scattering cross section. E † (r′, t′, n; r j ) is the adjoint energy field generated by the adjoint source F † (t, r″) which contains the information about misfit. E † is obtained as: since there is no energy pattern changed (Margerin et al., 2016). The adjoint source F † (t, r″) is derived from the misfit function as: where δ(r″ − r j ) is the Dirac function. Equation 9 shows how the match between model prediction and observed data enters the Fréchet derivative.
The iterative inversion starts with an initial model 0 = { 0 = 2 0 ( ), 0 = −1 0 ( ) } . The forward and adjoint fields are simulated with this initial model to calculate the scattering and absorption misfit kernels 0 ( ′ ) and 0 ( ′ ) based on Equations 6 and 7 which represent the gradients of the misfit function. To minimize the misfit function, the search direction { , } of the model is calculated from the gradients { , } using L-BFGS method (Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm). Note when k = 0, the steepest decent method is used instead of L-BFGS method. With the appropriate step length, the model m k of iteration k will be updated as m k+1 : where η is the step length. More details about L-BFGS and the choice of step length refer to T. Zhang and Sens-Schönfelder (2022). The final model is obtained by repeating the workflow described above until the misfit converges.

Experiments
To test the AET against real data we choose an acoustic experiment conducted in a reinforced concrete specimen at the German Federal Institute for Material Science and Testing (BAM). The sample has a size of 4 m × 5 m with a height of 0.8 m as shown in Figure 1 (Epple et al., 2020). In this paper, the orientation defined as shown in Figure 1 is for convenience to discuss and not the natural geographic coordinate. All directions in the following discussion refer to this definition of orientation. Nineteen ultrasonic transducers are embedded in the central layer of the specimen at 0.4 m height. The transducers serve as the energy sources of ultrasound with a center frequency of 60 kHz and as receivers. Both emission and recording of acoustic waves is laterally isotropic. This setup provides for 19 × 18 source-receiver combinations. The experiment has the following advantages for the present purpose: (a) due to the rather flat shape of the specimen and the placement of the transducers in its central plane we can restrict the energy propagation to the lateral directions and simplify the problem to 2D. (b) The boundary conditions of the lateral edges of the specimen can easily modeled using mirror sources. (c) The 10.1029/2022JB024972 6 of 21 embedded sensors that are located 30 cm away from the free surfaces reduces the excitation of surface waves which are not treated in our approach.

Data Processing
The ultrasound signals were recorded with sampling interval of 0.5 μs for a lapse time of 5 ms. Seven identical experiments were performed on three consecutive days in October between 7:00 and 8:00 a.m. An illustration of original data excited at source T0120 and recorded by receiver T0135 is shown in Figure 2a. The first 200 samples, that is, 0-0.1 ms precede the signal transmission and are recorded to control the noise level (Niederleithinger et al., 2018). The impulsive signal at 0.1 ms lapse time is visible on all sensors and is caused by electro-magnetic cross-talk between the cables of the high voltage source and the recording sensors. We use this signal to synchronize the pulse generator and the acquisition unit and determine the source time. Data is detrended by subtracting its mean (Figure 2b). The cross-talk is used to extract the envelope of the source signal and is then removed from the record (Figure 2c) which band-pass filtered between 60 and 120 kHz to remove the high-frequency noise ( Figure 2d). The envelope of the filtered signal is extracted using the Hilbert transform ( Figure 2e). Envelopes of the repeated experiments are averaged to obtain the final envelope ( Figure 2f) for the inversion. The same processing is applied to the cross-talk to obtain the final envelope of the source signal (Figure 2f inset).
From the processed envelopes we noticed that certain sensors systematically recorded smaller amplitudes than others, or excited less energetic waves. We attribute this to variable sensor coupling including the conversion between electrical and mechanical signals as well as the mechanical coupling between the transducer and the concrete. We estimate the coupling using the coda normalization method (Sato et al., 2012) which states that the signal envelopes in the late coda should be independent of location due to the equal distribution of elastic energy. We estimate one coupling coefficient for each transducer acting as source and receiver, separately by averaging the late coda envelope (3.5-4.7 ms) from the respective source or recorded at the respective station ( Figure 3a). Since the transducers act both as source and receiver, the coupling should have similar effects on both 7 of 21 the emission and recording. This is consistent with the observations in Figure 3a. The influence on the envelope data from the ith source to the jth station is eliminated by dividing by the corresponding values in Figure 3a. An illustration of the coupling effect is shown in Figure 3b. The blue and red curves indicate two combinations exchanging the source and the station, which should be identical due to reciprocity. However, the sensor coupling introduces a difference between two curves shown in Figure 3b but can be corrected using the coupling corrections ( Figure 3c).

Diffusion Model
As introduced in Section 2.2, the iterative inversion starts with an initial model. The density of the concrete is provided by Niederleithinger (2017)  We have no prior information about the scattering and absorption properties, for the concrete in the present experiment. A simple description of multiple-scattering and intrinsic attenuation of ultrasound in concrete is provided by the diffusion model (Anugonda et al., 2001;Ramamoorthy et al., 2004). The 2D diffusion equation describes the energy radiating isotropically from a source : The diffusion energy E D (r, t; r i ) at position r with the lapse time t is determined by the source energy E 0 , diffusion constant D and intrinsic factor Q −1 at the specific angular frequency ω. 1/4πDt and − 2 ∕(4 ) describe the shperical diffusion of energy in 2D and e −ωt/Q is the intrinsic attenuation term. r i is the position of the source while the distance between the source and receiver is = | − | . To account for the existence of boundaries that reflect the acoustic energy, E D (r, t; r i ) is summed for all mirror sources corresponding to r i (Y. Zhang et al., 2018).
Equation 13 consists of three terms in which lnE 0 is constant. To speed up the process, we separately estimate Q −1 from the later coda wave (3.5-4.7 ms) since the later coda wave is more sensitive to the intrinsic attenuation (T. Zhang & Sens-Schönfelder, 2022).
[ ln(4 ) + 2 ∕4 ] varies slowly in the late coda. Therefore, −ω/Q is easily estimated from slope of the logarithmic envelope in the late coda. Figure 5a shows the distribution of the estimated Q −1 values from all source-sensor combinations. The mean and median value of this distribution are both 0.003 that will be used to estimate D in the diffusion modeling and as initial model for inversion.
With the fixed value of Q −1 , the source energy E 0 can be extracted as the offset from the envelopes for each assumed diffusion constant D. An interval of [50, 5000] mm 2 ∕s with step-length 10 mm 2 /s is searched for the   Figure 2f to compare with the observable. For all source-receiver combinations, the distributions of D and E 0 are shown in Figures 5b and 5c. According to this distribution, we fix the source energy E 0 in this study to 12. The diffusion constant D does not directly correspond to the parameters used for the non-isotropic scattering in RTT. It corresponds to the transport scattering coefficient g* which is a version of g that is weighted by the cosine of the scattering angle θ. The relationship between D and g* is given as where * 0 is the average transport scattering coefficient that is defined as: g(θ) has been introduced as a function of scatter strength ɛ and correlation length a in Equation 2. Assuming that the correlation length a is uniform with a = 0.011 m (Anugonda et al., 2001) we calculate the values of ɛ corresponding to the estimated values of D using Equations 2, 14, and 15. The distribution of ɛ is shown in Figure 5d. We fix ɛ = 0.13 as background parameter describing the small-scale heterogeneity in the concrete specimen. Table 1 summarizes all background parameters estimated for the use with Equation 1.

Monte Carlo Simulation
The radiative transfer equation is solved using the Monte-Carlo method to simulate the energy propagating (T. Zhang et al., 2021). To account for the free surface boundary conditions in the Monte Carlo simulations the particles are reflected at the four sides of the model. In this study, 100 million particles are used for each simulation. The field generated by source T0120 in the initial model is illustrated in Figures 6a-6f. Although the algorithm allows us to simulate in models with spatially variable ɛ 2 (r) and Q −1 (r), here we only illustrate propagation in an uniform model with the background parameters given in Table 1. Note Figure 6 only shows the energy density ( , ) , while we actually simulate the specific energy density ( , , ) with information about the propagation direction.
The Monte Carlo method simulates a point-source in space and time. The simulation result is therefore convolved with the source wavelet and multiplied with the same source energy E 0 as diffusion model. Figure 6g shows a comparison of one observed envelope with the diffusion model and the MC simulation in the background model. The blue and red curves represent the energy simulated with the diffusion model and radiative transfer equation, respectively.

Imaging
Starting from the initial model with uniform parameters estimated with the diffusion approximation, we use AET to infer the spatial distribution of the strength of heterogeneity and attenuation. Both material properties influence the energy propagation causing a trade-off between changes in the scattering and absorption properties in a simul taneous inversion for both parameters as discussed in T. Zhang and Sens-Schönfelder (2022). For the ballistic wavefield, that is, the energy that propagates without being scattered, the effect of scattering and attenuation is identical -leading to the impossibility of discerning both effects with direct waves. But the trade-off also exists for arbitrary sub-segments of the propagation path of coda waves. Only the combination between the energy that propagates directly between two points in the medium and the energy that is scattered between these points allows us to resolve the trade-off since heterogeneity increases the scattered part of the wavefield at the cost of the direct part. This trade-off means that strong spatial differences of one parameter unavoidably map into the other parameter to some extent (Cormier & Sanborn, 2019;T. Zhang & Sens-Schönfelder, 2022).
However, the fact that the early coda is important to image the heterogeneity while the later coda is more sensitive to intrinsic attenuation (Calvet et al., 2013; T. Zhang & Sens-Schönfelder, 2022) helps us to separately invert ɛ 2 (r) 2π ⋅ 60 kHz 2.4 g ⋅ cm −3 4.475 m ⋅ ms −1 0.003 12 0.13 0.011 m

Table 1 All Background Parameters Estimated for the Radiative Transfer Equation
and Q −1 (r) using the early and later coda, respectively. In this experiment, we simply define the early and later coda intervals by 1.7-3.5 ms and 3.5-4.7 ms respectively as shown in Figure 6g and use these time windows to image the absorption and scattering structures successively.

Intrinsic Attenuation Inversion
We first focus on the intrinsic attenuation inversion with Q −1 (r) since the absorption influences the whole envelope. The later coda wave (3.5-4.7 ms) is chosen as the time window to evaluate the misfit function and the initial model −1 0 ( ) is uniform with −1 0 = 0.003 . The other parameters and the model of ɛ 2 (r) are all uniform based on Table 1 and remain constant during the inversion, meaning that only −1 ( ) is updated.
After 11 iterations of AET, the normalized misfit between the observed envelopes and synthetic data converges to 66% as shown in Figure 7a. The decrease of the misfit is very fast in the beginning since the initial model is uniform, slows down and stagnates from iteration 7. The benefit of iterative inversion as compared to a linear kernel-based inversion (Ogiso, 2019) is that the model is further improved after the first iteration based on the results of earlier iterations. The final inversion result is shown in Figure 7b. The distribution of Q −1 (r) shows a dominant first order structure with a maximum in the center and a symmetry in the west-east and north-south directions. The decrease toward the sides is not isotropic with the east-west direction showing faster decrease than the north-south direction. We will discuss the interpretation of this result in the next section.
The misfit is the integral of the differences between the observed and modeled results in the specified time window. However, we can also directly check the data fit of the envelopes. Figure 8 shows the data fit for Figure 6. (a-f) The snapshots of the simulated energy density field from T0120 (red star) at different lapse times. The scattering mean free path is 0.36 m and the mean free time is 0.08 ms. (g) The comparison among the envelopes recorded at T0135 (red inverted triangle) from the Monte Carlo simulation (red curve), the diffusion model (blue curve), and the real data from the concrete experiment (black curve). Note that the color scale range of each time in (a-f) is different and the energy density field has been multiplied with E 0 but not convolved with the source so that the values are not the same as the envelopes shown in (g). ZHANG ET AL.

10.1029/2022JB024972
11 of 21 some source-receiver combinations. The simulated envelopes in the final inverted model (red solid curves) are compared with the initial model (blue dashed curves) and the observed data (black solid curves). The locations of source and receiver in each combination are shown on the right side. For the north-south oriented combinations T0119-126, T0123-130, T0126-133, and T0130-137 which are located in the west and east, the envelopes of the inverted model become more similar to the observation compared with the initial model, as expected for a successful inversion that minimizes the misfit. This is caused by the decrease of Q −1 along the western and eastern sides of the model. There are no significant improvements for station combinations T0120-122 and T0134-136 because already the initial model fits the observations reasonably well in these areas and the model update during the inversion is marginal. Envelope fits of the combinations T0124-125 and T0131-132 that transect  12 of 21 through the whole specimen do not improve. In fact the fit of these long distance east-west combinations slightly degrades in favor of significant improvements of other pairs.

Scattering Inversion
Although the early coda wave is more sensitive to scattering, scattering inversion can benefit from using a more reasonable model of Q −1 to suppress the influence of the absorption. In this step, we employ the inversion result shown in Figure 7b as the model of Q −1 (r) and keep it constant throughout the inversions for ɛ(r). The initial model of ɛ(r) is uniform and we use the earlier time window with lapse times 1.7-3.5 ms (cf. Figure 6).
Nine iterations were conducted until the normalized misfit converged to 77% which is shown in Figure 9a. Note that although the simulation in the initial model in Figure 9a is the same as the last one in Figure 7a, the absolute value of misfit is not since the time windows are different. Figure 9b shows the inversion result of ɛ(r). The dominant value of it is about 0.14 which is a little higher than the initial uniform model 0.13. The inferred distribution of heterogeneity has a more complex structure than the attenuation structure. Stronger scattering is inferred in two areas at the western and eastern boundaries and also in one anomaly of higher value in the south at about y = 2 m. An elongated features extends from the northern to the southern edge at about y = 3.4 m. A very low-value anomaly that indicates reduced heterogeneity is located in the north-east corner. Interpretations are discussed in the next section.
The data fits are shown in Figure 10. Similar to what we discussed in Figure 8, the different-distance combinations are compared in the early coda waves time window. The inversion result is dominated by the short-distance combinations which achieve a significantly improved data fit during the inversion. The medium-and long-distance combinations do not improve clearly. Note that the y-scale of the graphs in Figure 10 is variable and combinations T0124-125 and T0131-132 have far smaller amplitudes.

Misfit Evolution
In Section 4, we have described two successive inversion runs for ɛ 2 (r) and Q −1 (r) using the early and later coda, respectively. We start with uniform models of both parameters, first update the model of Q −1 (r) only, and then fix the Q −1 (r)-model and continue to update ɛ 2 (r). The time windows of the misfit are chosen to use only the later coda for intrinsic attenuation inversion and only the early coda for scattering inversion because of their sensitivities. Of course, the misfits of both time windows varied in both inversions. Figure 11 shows the evolution of the misfits of both time windows for the whole inversion process. The red and blue curves indicate the misfits of the later and early coda, respectively. The solid parts of the curves show the misfits that are optimized for during The whole inversion is separated into two periods shown in Figure 11. In the first period when we only update Q −1 (r) (red domain), the misfit of the later coda (the red solid curve) decreases since the misfit kernel is based on this time window. Reasonably, with the improvement of the attenuation model the misfit of the early coda (the blue dashed curve) decreases as well although it is not used to guide the inversion. During the subsequent updating of ɛ 2 (r) (blue domain) the misfit of early coda time window continues to decrease since it is used to calculate the adjoint source. On the contrary the misfit of the late time window which is not used in this step re-increases Figure 10. The data fitting of different combinations (illustrated on the right side, red star and blue inverted triangle are source and receiver respectively) in the early coda wave (1.7-3.5 ms). The blue dashed and red solid curves indicate the envelopes simulated in the initial model and the inverted model shown in Figure 9b, respectively. The black curve is the real data from the concrete experiment. slightly which is not surprising since this time window was already optimized for in the Q −1 (r)-inversion and does not inform the ɛ 2 (r)-inversion. However, the misfit change in the second run is dominated by the decrease of the misfit in the early time window. Using both time windows together to guide the second part of the inversion run would possibly have damped the misfit increase in the late time window, at the expense of smaller improvements in the early time window.

Interpretation
We begin the discussion with an interpretation of the inferred attenuation. The attenuation anomaly (Figure 7b) is symmetric with respect to west-east and north-south axis in the center of the specimen and appears to be affected by some large scale influence on the specimen rather than internal small scale differences. Three processes could globally affect the specimen and result in a perturbation with the symmetry observed in the attenuation structure: (a) diffusion of humidity, (b) temperature changes and (c) stress distribution.
To investigate this hypothesis we make use of supplemental instrumentation. Additional to the 19 ultrasonic sensors, there were four temperature sensors embedded in the concrete specimen (shown in Figure 1) which measured the internal temperatures on three consecutive days in the morning between 6. and 8 a.m. as shown in Figure 12a. This experiment was conducted during a phase of decreasing temperatures in autumn. The temperature at each sensor decreased during the three successive days but the sensors maintained rather constant offsets from one another. The central sensor T0128 shows highest temperatures compared to the sensors closer to the rim. Smallest temperatures are observed in the corner of the specimen at sensor T0137 while intermediate temperatures are observed along the sides. We use the temperature measured on 28 October at these four sensors to obtain an idea of the temperature distribution within the specimen. We therefore use the geometric symmetry of the sensor locations to interpolate the observations throughout the whole concrete in 2D using adjustable tension continuous curvature splines by Generic Mapping Tools (Smith & Wessel, 1990;Wessel et al., 2013). The resulting temperature distribution within the concrete is shown in Figure 12b. This is clearly a rough estimate of the internal temperature distribution, but it shares clear similarity with the inversion result of Q −1 (r) shown in Figure 7b.
It has been demonstrated that the temperature changes of the concrete can result in the velocity perturbation but the sensitivity is only about 0.05-0.15%K −1 (Epple et al., 2020;Larose et al., 2006;Niederleithinger & Wunderlich, 2013). Since the maximum temperature change during the experiment is only 0.5 K(°C) the observed temperature changes will thus have a negligible influence on the propagation velocity and thus leave the envelopes unaffected which warrants the assumption of uniform and constant velocity in this experiment. We did not find conclusive evidence in the literature for the influence of temperature on attenuation in concrete or similar aggregates (at the present temperature (Zong et al., 2020)). The influence of humidity on attenuation has been clearly documented by a number of authors (Clark et al., 1980;Green et al., 1993;Tisato & Quintal, 2014). Unfortunately in-situ observations of humidity are not available to us and the specimen is insulated from the sides and covered for protection against rain so that the humidity might be more or less uniform in the volume.
A distribution of absorption with a very similar symmetry pattern in a sample of comparable size was found by Liu and Guo (2005). These authors imaged the attenuation in a reinforced concrete block under the highway bridge pier cap which had a size of 6 m × 8 m with a height of 1.5 m. Using direct waves Liu and Guo (2005) inferred an inverse intrinsic quality factor of 0.0063 in the center of the block. This value is close to our result 0.0045. Toward the sides of their block, attenuation increases seven times while it decreases four times in our results.
Different from absorption, the heterogeneity of the medium appears to be governed by internal structure rather than an external influence since the inferred distribution is much more structured. Figure 13 shows the construction drawing of this concrete specimen. The strongest anomaly of increased heterogeneity is found at the western edge of the specimen. This area corresponds to a volume of the specimen that was cast with a different kind of concrete (salt concrete: 1,600 × 1000 × 250 mm). Here, salt was added to the concrete mix to be able to provoke rapid corrosion of rebar at a later stage. As the concrete was poured separately by a different team and cured under different conditions, a different density and porosity can be expected. We interpret the increased scattering inferred in this region to be caused by the different properties of the salt-concrete.
The second prominent area of increased heterogeneity located in the east does not directly correspond to model features from the construction plans. During the installation of the embedded sensors an anomaly was detected in this area. While the calculated quantity of grout was sufficient to completely fill the boreholes in all other locations, almost three times the amount was required for refilling the borehole of sensor T0132. It can therefore be assumed that cavities were unintentionally created in this area during concreting, which now contribute to the increased scattering.
Before this experiment, there were three heating cartridges inserted in the south, east and northeast (Heating Cartridge A, B, and C respectively in Figure 13). Heating Cartridge A had been used to heat the concrete to 510°C (Niederleithinger, 2017) while the other two had not been activated. The concrete after high-temperature heating generated thermal cracking and stress changes (Hager, 2013) that increase scattering.
Three autoclaved aerated cube concretes with size of 0.3 m, four horizontal plastic pipes and one vertical clamping channel are also embedded in the concrete. Structures of these size are out of the inversion resolution but can also affect the scattering to some degree.
A very prominent anomaly that is left to be discussed is the low-δɛ anomaly close to sensor T0122. This anomaly is located right at the boundary of the inversion domain and converges toward extremely low values of heterogeneity, that is, locally homogeneous material. Its location directly on the boundary close to a corner of the model leads us to the interpretation as an artifact. Fitting envelopes of waveforms always requires significant averaging. In theory this averaging should be achieved by repeated observations in statistically identical realizations of the experiment. In reality there is only a single specimen and the averaging is realized on the one hand based on ergodicity by using long time windows for the comparison between observations and synthetics and on the other hand by using multiple source and receiver combinations. While the effect of long time windows is the same everywhere in the sample the averaging by different sensor combinations is not. The reflecting boundary conditions reduce the effective averaging by a factor of two along the edges and by a factor of four in the corners. A prominent wiggle in the waveform that can coincidentally originate from the constructive interference of scattered waves results in a strong pulse in the envelope (cf. Figures 6, 8, and 10). Such a pulse can push the inversion into a certain direction and cannot effectively be compensated by other sensor combinations with sensitivity to the same location since the mirror sources have identical waveforms.

Resolution Test
Different tools exist to study the capabilities of the combination of a measurement setup and an inversion method. Checkerboard tests (Lévěque et al., 1993) use a periodic pattern of variable wavelengths to infer the minimum size of a feature to be resolved in different parts of the domain. Analytical approaches use the sensitivity of the misfit function to changing perturbations (the Hessian) at the different locations in the domain to estimate the resolution capabilities (Fichtner & Trampert, 2011).
We take a different approach for the following reasons. Since we use reflecting boundary conditions in a domain with a regular distribution of sensors we can assume that also the resolution capabilities are rather uniform which would limit the value of a checkerboard test. The analytic approach using the Hessian is either computationally very expensive or requires further development, that is beyond the present scope. Here we ask the question: What would the inversion obtain if the structures were as we interpret it from the actual imaging. Technically this question is answered by inverting a simplified version of the obtained result that contains all structures which are regarded as relevant and interpreted. This approach is often used in tomography to confirm that the interpreted structure could indeed be resolved by the imaging (Jiang et al., 2014;Koulakov et al., 2009). For nonlinear problems such statements are more useful than theoretical values of resolution length in a homogeneous background model.
In our resolution test we set up a test model, generate synthetic data by forward modeling and investigate how well the inversion procedure is able to recover the original model. The forward simulation that computes the synthetic envelopes in the test model uses the same Monte-Carlo method that is also used in the inversion. The inversion of the synthetic data is then performed using exactly the same procedure and parameters (sensor location, time windows, starting models, etc.) as in the original inversion of the real data in Section 4.
The test model of ɛ(r) is based on the construction plans of the concrete specimen and the inversion result. The background value of ɛ(r) is designed not to be the same as the initial model but taken from the inversion result as 0.14. Figure 14a shows the input model for the resolution test that contains the structures obtained in the inversion 17 of 21 and some small elongated anomalies along the locations of channels and reinforcement bars in the specimen. The input model of Q −1 (r) (Figure 14b) is based on the temperature distribution shown in Figure 12b. Panels (c) and (d) of Figure 14 show the resulting outputs of the synthetic inversion test. The output of ɛ(r) shows that the background value is recovered well although it was different from the initial model. The three larger anomalies are localized well, but their shapes are not recovered in detail due to limitations imposed by the number and the setup of sources and receivers and the intrinsic smoothing of imaging with the envelope information, only. Likewise the thin elongated anomalies are not resolved as could be expected from the locations of the 19 sensors of which only three are not arranged along the rim of the specimen. The incorrectly inferred shape of the anomalies is connected to their peak amplitudes which are partially overestimated during the inversion. Since the scattered energy depends to first order on an integral scattering strength of the anomaly higher values in the centers of the larger anomalies compensate for the lower strength along the edges of these anomalies. The inversion result of Q −1 (r) recovered the input structure well. However, decay in the north-south direction is underestimated and the peak anomaly is overestimated.
From this test, we conclude that the first order features interpreted from the imaged attenuation and scattering structures would indeed show up as observed in the results. Due to ambiguity and limited resolution we cannot exclude that smaller anomalies are present in the specimen.

Solid Earth Applications
AET is applicable to elastic waves in the multiple scattering regime. The existence of this regime depends on the characteristics of the heterogeneity and on the wave length of the elastic waves but it can be found in ultrasound laboratory experiments with concrete and in earthquake records at high frequencies. The frequency in our experiment is 60 kHz corresponding typical frequencies of about 1 Hz in seismological investigations of the scattered 18 of 21 wavefield. Since the acoustic velocity in concrete does not significantly differ typical seismic velocities in crustal rocks, the size of the block of concrete would correspond to a field scale of about 300 km which is rather typical for studies of wave scattering at volcanoes (De Siena et al., 2013Del Pezzo et al., 2016).
The setup of the sensors shown in Figure 1 results in source-receiver distances ranging from 0.85 to 5.56 m which translates into reasonable epicentral distances up to a few hundred kilometers. Station numbers and spacing can be expected to be superior in field applications where several tens of stations can often be found used within an area of a few hundred kilometers radius. The degree of heterogeneity in concrete, however, is comparably high as can be readily recognized from the waveform shown in Figure 6. A direct wave is not visible in this waveform due to scattering and confirmed that the situation in concrete corresponds to those parts of the Earth that exhibit strong heterogeneity such as fault zones or volcanoes.
Besides the length scale, our assumptions of 2D geometry and acoustic wave propagation may have implications for solid Earth applications. The 2D geometry is similar for investigations of scattering in the crust due to the more homogeneous mantle (Sens-Schönfelder et al., 2009). In terms of computational efficiency, the 2D case is similar to 3D scattering in spherically symmetric settings like problems. If heterogeneity varies only in a single dimension, for example, with depth, the computational demands are even smaller. Lastly the assumption of acoustic scattering is often made also in field application to the Earth because scattering leads to a stable balance between P and S energy in the scattered wavefield called equipartition (Hennino et al., 2001). In this state the wavefield is dominated by shear waves and scattering can be treated as single mode waves in the acoustic approximation. Another interesting application that allows to neglect mode conversion in the scattering process is the investigation of early teleseismic phases such as P diff coda (Earle & Shearer, 2001) and the PKP precursor Thomas et al., 2000). For these wavefield the large travel time filters out the slower S-wave energy. A very promising application of AET for the future is the construction of a 1D Earth model for strength of heterogeneity in analogy to existing reference models for velocity (Dziewonski & Anderson, 1981) or attenuation (Montagner & Kennett, 1996). Even improved starting models exist for the heterogeneity that have been derived with less inversion strategies (Bentham et al., 2017).

Conclusions
This research presents the analysis of an acoustic experiment conducted in a 4 m by 5 m large concrete specimen equipped with embedded acoustic sensors. We applied AET to image the distribution of small-scale heterogeneity and intrinsic attenuation inside the specimen. To interrogate the structure below the resolution limit of conventional tomography, AET was proposed to invert for the statistical properties of the small-scale heterogeneity as complementary information to the deterministic structures that can only be imaged at larger scales. Although AET had been successfully tested in numerical experiments, the application to experimental data in the present paper increases confidence in the methodology in view of further applications to seismic imaging of the Earth.
We performed this experiment with ultrasonic transmission from embedded transducers in reinforced concrete in analogy to seismic wave propagation in the Earth. The data recorded by 19 transducers are compared with simulations of energy propagation based on the radiative transfer equation. This forward problem is solved by modeling the 2-D multiple nonisotropic scattering in an acoustic medium with spatially variable heterogeneity and attenuation using the Monte-Carlo method. The misfit between the observed and modeled envelopes is minimized by iteratively updating the model with the adjoint method. The whole workflow of AET for the real data is introduced including the processing of the data and the investigation of background values with the diffusion model. The fluctuation strength ɛ and intrinsic quality factor Q −1 respectively representing the spatial variability of scattering and absorption are separately inverted from different time windows. On the one side, the absorption inversion result shows a strong point-symmetric geometry which we interpret as some large-scale spatially variable in the specimen, but without a direct evidences for the causative process, for example, temperature, humidity or stress.
The inverted distribution of scattering properties shows a more complex structure that can-to some extent-be interpreted in terms of the known internal structure of the test specimen. The largest anomaly of increased heterogeneity corresponds to a volume containing salt-concrete. Other anomalies are not as clearly linked to the known features of the concrete and a strong anomaly of decreased heterogeneity exists at the edge of the specimen that is interpreted as an artifact from envelope fluctuation that are insufficiently averaged at the reflecting boundaries of the model domain.
Despite obvious room for improvement in terms of spatial resolution and power to resolve the trade-off between scattering and attenuation the present results are encouraging. The spatial variability of attenuation and scattering strength improved the data fit by about 35% when averaged over both time windows. This number appears small but cannot directly be compared to improvements known from waveform inversion. Two effects contribute to the limitation of the data fit. First the observed envelopes are obtained in a real experiment and cannot be averaged over an ensemble of test specimens and thus show fluctuations introduced by the interference of scattered waves that cannot be fit. Second also the simulated envelopes contain additional fluctuations from the Monte-Carlo type simulation.
Future investigations to test the performance of the AET on real data will have to include dedicated test specimens with known scattering and attenuation properties. Even though the present concrete block with the embedded sensors was well suited for an application of AET, it was already cast and the different types of concrete could not be analyzed separately to obtain ground truth. An important field of application for the presented approach is the monitoring of medium perturbations with coda waves (Sens-Schönfelder & Brenguier, 2019). The spatial sensitivity of coda wave based monitoring depends on the distribution of heterogeneity (Kanu & Snieder, 2015) and can thus be improved with the presented method. We hope that AET will contribute to non-destructive testing of civil engineering structures and investigations of wave propagation in the Earth.

Data Availability Statement
The code for Monte Carlo simulation, scripts used to calculate the misfit kernels and the processed data of the laboratory experiment can be accessed at https://doi.org/10.5281/zenodo.7152278.