A Cyclic Multi-Stage Implementation of the Full-Waveform Inversion for the Identiﬁcation of Anomalies in Dams

: For the safe and efﬁcient operation of dams, frequent monitoring and maintenance are required. These are usually expensive, time consuming, and cumbersome. To alleviate these issues, we propose applying a wave-based scheme for the location and quantiﬁcation of damages in dams. To obtain high-resolution “interpretable” images of the damaged regions, we drew inspiration from non-linear full-multigrid methods for inverse problems and applied a new cyclic multi-stage full-waveform inversion (FWI) scheme. Our approach is less susceptible to the stability issues faced by the standard FWI scheme when dealing with ill-posed problems. In this paper, we ﬁrst selected an optimal acquisition setup and then applied synthetic data to demonstrate the capability of our approach in identifying a series of anomalies in dams by a mixture of reﬂection and transmission tomography. The results had sufﬁcient robustness, showing the prospects of application in the ﬁeld of non-destructive testing of dams.


Introduction
Tomography is a well-established technique in the non-destructive testing of bodies and structures. It is based on a technique in which signals of a particular physical domain are sent through a system. Computer programs evaluate the measurements of the reflected and transmitted signals and visualize the interior conditions, e.g., by the prognosis of the material distribution inside. In the analysis of structures in civil engineering, most tomographic methods are based on mechanical waves. This technique is well known, e.g., in the field of geotechnical exploration via seismic tomography [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. The wave equation is the basis, and by an inversion, the layers of different densities and stiffnesses are detected by reflection tomography. The current paper puts forward the first steps in the adoption of this technique for the assessment of the structural condition of dams. As is well known, dams are important structures employed for the generation of hydro-electricity, the provision of water supply, and flood defense. However, a large number of the dams that are in operation today were built some decades ago, and their structural condition is not properly known. The effects of damage, aging, chemical reactions, and internal erosion may have deteriorated the material properties. It is of the utmost importance to assess the structural condition carefully to make further trustworthy predictions of the structure's reliability [15]. Unfortunately, the conventional inspection methods are generally time consuming, expensive, and not precise enough; see, for example [16]. To tackle these challenges, a series of methods was proposed by [17][18][19][20][21][22][23][24][25][26][27][28]. These methods were designed to find anomalies in technical structures, mainly for coupled hydro-mechanical systems or strongly related systems.
Using classic ultrasonic exploration methods, such as the synthetic aperture focusing technique (SAFT) or the total focusing method (TFM), will not provide the desired penetra-tion depth, and only areas close to the surface might be assessed, as all the information of the ultrasonic waves is not used. However, with the recent technological improvements in computational capabilities and transducer technology, the usage of the FWI originating from seismic wave probing may be extended to structures of smaller scales than those typical in sub-surface exploration. The challenge in dealing with this is that visibly shorter wavelengths occur, which, first of all, require efficient solvers of the wave equation itself. Moreover, this is a critical effect using higher excitation frequencies that results in an increase in the ill-posedness of the parameter identification problem (finding anomalies in dams). A straight-forward application of FWI seems, therefore, impossible, and algorithmic extensions need to be made, which is the main contribution of this work.
FWI, initially proposed by [29,30] and further developed in the last few decades, can be regarded as a modern seismic imaging technique. As the subsurface properties do influence waves, in particular, pressure and shear waves, they can be utilized to probe heterogeneous domains. FWI is actually an inverse parameter identification problem based on the wave equation. It may consider recorded signals as the input for a fixed time according to some dynamic excitation. The information in the signals includes the effects of reflection, refraction, diffusion, and amplitude. Thus, they are comparatively rich in information to infer the mediums' material parameters. The approach of this method, based on the minimization of the misfit between recorded and numerically generated data (or a regularized version of them), allows for the estimation of the parameters influencing the propagation of the seismic body waves' velocities and density [31]. From the wave velocities, the Lamé parameters, defining the stress tensor, can be derived. FWI offers unique advantages in terms of flexibility and increased resolution (e.g., compared to linearized approaches) and is, hence, capable of imaging the arbitrarily heterogeneous compression and shear wave velocity profiles of objects of any shape [2].
Based on the successful application of full-waveform inversion (FWI) in geotechnical exploration and non-destructive testing (NDT) [32][33][34][35], especially with subsurface and engineering structures, the presented method extends this with an application to dams. In addition to this, the algorithm was formulated as a cyclic multi-frequency stage inversion exploiting the potentials of FWI in a robust extraction of the complete waveform information to produce high-resolution images.
Of course, inverting wave equations is itself computationally demanding. Embedding FWI into a cyclic multi-frequency-stage algorithm increases the demands even more. On the other hand, no other known inspection tool offers such a high accuracy that the increased computational demands are justified; compare, e.g., [16,36].
In the upcoming sections of this paper, we take an in-depth look at the following: • Materials and methods: the formulation/algorithm of the forward model and its boundary conditions, the inverse model, and the inversion process. • Numerical simulation: the geometry of the dam, the acquisition geometry, the wave parameters, and the cyclic multi-stage implementation under different conditions. • Results: FWI results for different acquisition geometries, cyclic multi-stage inversion, and for different levels of noise. • Conclusion: a summary of the applicability of our approach to identifying damages in dams.

Materials and Methods
The paper at hand proposes the characterization of a dam's structural properties by applying FWI. Therefore, both the forward and inverse model are presented in the sequel, including some computational details and the extension of the existing theory by introducing a cyclic multi-frequency-stage inversion scheme.

Forward Model
The wave equation belongs to the class of hyperbolic partial differential equations. The equation is merely describing the propagation of a wave through a medium. Its solution gives the wavefield at given time and space. The full-waveform inversion is mainly making use of the wave equation formulated as the acoustic wave equation or the elastic wave equation.
The system of partial differential equations in Equation (1) is characterizing the elastic waves propagating in a 2D isotropic linear-elastic medium. Furthermore, only non-zero displacements in the horizontal direction and depth (x-z-plane) are described [1,37].
where i and j represent the directions x or z and perfectly matched layers are used as the boundary conditions. v x , v z are the particles velocities, t is the time, λ, µ are the Lamé parameters, ρ is the density, σ xx , σ xz , σ zz are the stress tensor components and f x , f z are the directed body force. Accordingly, the particle velocities can be calculated from the Lamé parameters using The elastic wave equation is solved numerically for the particle velocities v x , v z , stresses σ xx , σ xz , σ zz , density ρ and Lamé parameters (λ, µ) using a time domain 2D finite difference scheme. It is discretized in time and space, with each of the parameters placed on a staggered grid [38,39]. The discretization needs to satisfy certain spatial and temporal sampling conditions in order to avoid numerical artifacts and instabilities in the finite difference calculation.
The grid spacing dh is chosen by the number of grid points for the smallest wavelength λ w min in the system. This approach is undertaken to prevent grid dispersion in the simulation. Furthermore, the Courant-Friedrichs-Lewy criterion (CFL) [40] needs to be fulfilled in order to stabilize the simulations, and the time stepping dt needs to take less than the time needed for the wave to travel between two neighboring grid points with distance dh. Due to the complex free surface of the dam, we apply the Improved Vacuum Formulation (IVF), proposed by [41] at the air-medium boundary (i.e., top) and the convolutional perfectly matched layers (C-PML) are applied at the other boundaries to ensure stability and accuracy of the results.

Inverse Analysis
The main part of the inverse analysis is the definition of the forward operator, also called the parameter-to-solution mapping.
In this publication, it uses the following notations Here, X is a finite dimensional space of parameters (describing material distributions, λ(x, y), µ(x, y) and ρ(x, y)), Y is a finite dimensional space of time-dependent measurements. Furthermore, m and u are denoting the model parameters (material distributions) and the model responses, i.e., wavefields at the sensor locations.
The main concept behind the inverse analysis is formulated through the solution of the following nonlinear equation: where δ comprises measurement errors and model uncertainties. As the solution of Equation (4) may not exist as the measured data u δ are not in the range of F, we follow a least-square solution approach, as given in Equation (5). As a measure for the misfit between modeled and field data, an objective function C f (m) based on the L 2 -norm of the data residuals is therefore defined: where the excitation of the waves is achieved with frequency ω l . We use an iterative procedure which results in a discrete gradient-type procedure where m n are the iterates at iteration n, µ n is the step length, is the gradient and H n an approximation of the Hessian matrix at iteration n. An inexact parabolic line search is employed when estimating the µ n [42], while the time-domain gradients are calculated using the adjoint-state method as presented in [43][44][45][46][47]. For the Hessian, either the identity is chosen, which renders the algorithm a steepest-descent method, or we choose an approximation according to a Broyden-Fletcher-Goldfarb-Shanno (BFGS) updating scheme.
The gradient direction is computed by solving the stress-displacement (σ − u) elastic wave formulation adjoint problem in Equation (7). Here, the perturbation caused by anomaly in the data space is back propagated from the receivers into the unperturbed medium [48] ρ The initial conditions are u i (x, t = 0) = 0, ∂u i ∂t (x, t = 0) = 0 and the final conditions are δσ ij (x, t = T) = 0 and δ δx j δσ ij (x, t = T) = 0.
The full-waveform inversion procedure is implemented in the DENISE (subwavelength DEtail resolving Nonlinear Iterative SEismic inversion) 2D time domain C-based code [49] and shown by the flow chart in Figure 1.

Calculation of residuals
Cyclic-multi-stage algorithm with increase and decrease in the excitation frequency (see Figure 2 Inherent to all inverse problems in engineering are the "ill-posedness" and the nonlinearity of the problem, which need to be attenuated. In the case of FWI, the first step is to formulate a good starting model, where the model properties are close to a global minimum. Therefore, the parts of the long wavelengths from the "true" model should be represented by the starting model as satisfactorily as possible. Additionally, pre-conditioning or regularization should be applied. A reliable tool was found when carrying out the inversion at different stages, i.e., for different excitation frequencies, to cope with the non-linearity and ill-posedness of the problem. Each stage is carried out with a certain frequency, and according to [50], the low frequencies are more linearly related to the model perturbations compared to high frequencies. Thus, the inversion stages start with lower frequencies and increase to higher-frequency components of the data. The inversion output of each stage is used as the input for the next stage; thus, regions of high and low particle velocities are being captured. As mentioned, higher stages have the advantage of providing a better resolution; however, stability of the inversion deteriorates even with careful and early stopping. We observe divergence at many positions in the domain, which finally would lead to incorrect interpretations. To mitigate this and to further regularize, the process is repeated by going back to lower frequencies and then increasing again. This results in a cyclic multi-stage algorithm, which embraces the updating of the model parameters by setting borders in the frequency range; described here as cycles. A graphical sketch of the idea of this algorithm is given in Figure 2. As shown, on each intersection with a horizontal line, a FWI is run with a different excitation frequency. Frequencies increase and then decrease in one cycle. Repeating this provides chances to obtain high-resolution identification and keeps the effects of noise and resulting artifacts in the images low. Thus, the initially strongly ill-posed problem could be regularized sufficiently so that it might be applied in practice. In particular, for the first cycle(s), it is recommended not to probe with the highest excitation frequencies. The last cycle should be followed by an increase towards the maximal frequency.

Numerical Simulation Example
For the numerical simulation, a dam model (true model) is used, and the material properties are defined as seismic velocities and densities that are summarized in Table 1 and the true, the initial and the residual model are depicted in Figure 3. The distribution of the model parameters is given here, where the defects which ought to be detected are not contained in the initial model. The residual model is giving the difference between the true and the initial model.  The wave source and sensor/receiver setups are shown in Figure 4. They are used to collect the required waveform data to perform the inversion. Hence, changes in source wave properties due to reflection, refraction, diffraction, etc., are detected by the sensors. Three different acquisition setups are considered ( Figure 4) for the FWI.
Setup 1 consists of 28 excitation sources 65 m away from the dam surface and a total of 165 sensors/receivers on the dam and the reservoir bed/foundation. Receivers are placed at an interval of 0.5 m on the dam and 1 m on the foundation. Setup 2 is similar to setup 1, in the way that the same number of excitation sources is employed; however, in total, 135 sensors are used. The receivers on the dam foundation are placed within 35 m of the dam's face. Lastly, the same numbers of sources and receivers in setup 2 are applied to setup 3; the difference here is the distance of the excitation sources from the dam's face. The excitation sources are closer to the dam's face (i.e., 35 m away) than in the other setups.
The acquisition geometry can be set up in practice by placing ocean bottom cables (OBC) on the reservoir floor using boats or remotely operated vehicles (ROV). These four-component receivers consist of 3-axial geophones to measure the wave speeds and a hydrophone for the pressure wavefield. In addition to these, an array of geophones can also be installed on the downstream part of the dam. Similarly, the wave sources are placed by boats/ROV and can be changed after every shot or to adapt to different acquisition geometries. The OBCs are preferred for this setup, since the effect of noise from surrounding activities or water waves are minimized this way. The true model of the dam contains anomalies both in the structure of the dam and in the foundation. The anomalies represent weaknesses in different regions and of different sizes. In the dam structure, the top anomaly represents a 20% reduction in the parameter values V p , V s and ρ, the second weakening has a 10% reduction and the third one (i.e., crack zone) possesses a 30% reduction. All the anomalies in the foundation are reduced by 70% in the material parameter values. The entire simulation domain has a height of 48 m and a length of 128 m. This includes the dam structure with 28 m height, the crest of the dam is 4.5 m wide and it measures 19 m at the bottom. The water reservoir has a size of 30 m times 80 m, the foundation is modeled with 10 m of depth and the surrounding air is also considered, but there is no change in the air parameters included in the inversion process (they are fixed). A 2D finite difference scheme on a staggered grid is used to solve the elastic wave equation. The computational domain is replaced by a 2.5 m (equal to 10 grid points) thick border with convolutional PML (C-PML) on both sides and at the bottom of the domain. Owing to the high number of the systems' global degrees of freedom (i.e., identifying three parameters for each node in a domain with 512 × 192 = 98,304 nodes), the only meaningful strategy for retrieving a solution is via an adjoint state method embedded into a Quasi-Newton method. Here, the quasi-Newton limited memory Broyden-Fletcher-Goldfarb-Shanno (l-BFGS) method is used [42,51,52].
Since the FWI, based on local optimization, is sensitive to the initial model, an initial model containing the long-wave part of the model to be resolved is required to enable convergence to the global minimum. In contrast with geophysical exploration, which requires a series of field experiments to obtain information for the initial model, a priori information on the parameters of the dam material can often be obtained, since most of these structures are man-made. Therefore, the material properties can be used as the upper limit during construction, as in this case. A closer look at Figure 3 shows that the true model in Figure 3a and the initial model in Figure 3b differ only in the inclusion of regions with deteriorated material properties in the true model (see Figure 3c). Therefore, the "as-built" material properties of the dam are taken as the initial model. To study the dam material, spike wavelets are emitted successively from each of the sources (e.g., airguns). A series of low-pass filtered spikes are used here, as shown in Figure 4. For each excitation, the source wavelet as shown in Figure 5 propagates for a duration of 0.1 s, and the time interval (dt) used for the analysis is 3 × 10 −5 s, which satisfies the Courant-Friedrichs-Lewy criterion for stability. The signals recorded by the 2D receivers as a result of the wave propagating in the medium are shown in a Figure 6. This seismogram corresponds to the acquisition setup 1 (cf. Figure 4a). Similar seismograms can be recorded for setups 2 and 3, where the number of traces is reduced, because the number of sensors is decreased. These synthetically generated seismograms are used as field data for the inverse analysis (i.e., FWI).

Results
The FWI is carried out in multiple frequency stages, starting from lower frequencies to higher frequencies (i.e., 0.4 kHz to 1 kHz) to alleviate the non-linearity and the ill-posedness of the problem. The result of the previous inversion with lower frequency is used as the starting model for the following inversion, where a higher frequency is applied. Therefore, with an increasing number of inversions (i.e., increase in frequency per stage), the location and boundaries between different particle velocity regions are delineating, and thus, the quality of the inversion output is enhanced.
Additionally, the cyclic multi-frequency approach introduced in Section 4.3 is improving the inversion results as the cyclic increase and decrease in the frequency band dampens unwanted artifacts.

Optimal Acquisition Geometry Selection
The different acquisition geometries used are shown in Figure 4. The effects of the acquisition setups on the inversion quality shall be investigated and, therefore, the inversion results for the final stage of each setup are shown in Figure 7 for each of the model parameters. Additionally, the residuals between the inversion results and the true model are given in Figure 8. Overall, the damage in the dam structure can be easily identified on visual inspection. Anomalies in the foundation of the dam were not correctly identified. The best way to identify the anomalies/damage in the dam structure and foundation is to use the parameter V s instead of using one of the other two material parameters V p or ρ. The comparison of the cost function value for each of the setups after 40 iteration steps in Figure 9 shows that acquisition setup 3 provides the best quality on the inversion results. One inversion cycle stops and proceeds to the next inversion stage when the relative misfit change between the last two iterations is below a predefined value, which is set to 0.01 in this case. As seen in Figure 9, the transition to the new frequency level may cause an increase in the cost function value at the first iterations of the new stage. The superiority of setup 3 is also confirmed by examining the cross-sectional profile cuts through the regions of damage/anomalies in the dam (see Figure 10). Damage areas saturated with water can be identified by superposition of the models V p and V s . Such zones are made visible by looking at the properties of V p , which are capable of traveling through water (at slower speeds) and Vs which do not travel through fluids but produce better resolution for damaged zones.

Influence of Noise (Disturbance in Obtained Data) on the Reconstruction Quality/Error
To test the robustness of the proposed FWI method for dams, the effects of random noise on the (synthetically generated) measurement data are studied. This investigation is also necessary, as it tends to factor in some of the uncertainties which may result from the model quality, equipment reliability in data recording, data handling and other factors on the field, such as vibrations resulting from the dams operation and other unknown sources. These can be captured in practice by calibrating/tuning the sensors/receivers (especially on the downstream of the dam) to these "ambient" vibrations and filtering them out during analysis. The synthetically generated field data are corrupted by different levels of Gaussian noise before carrying out the inversion. For each of the sensors used in the acquisition setup, the recorded data are corrupted by 1%, 2%, 5% and 10% noise. A comparison of the cost function after each iteration for different noise levels, and also a comparison of the final cost function values for the noise levels are obtained in Figure 11. As for the comparison of the different setups, the relative misfit change is set to 0.01 before the inversion proceeds to the next inversion stage or is aborted for the last and final stage. Observations in Figure 11a show that the cost function trend for higher noise levels is shorter compared to lower noise levels. This behavior results from noise-induced velocity values, which are above the threshold set to ensure the CFL stability criterion; hence, the early convergence to a local minimum. Furthermore, in Figure 11b, the cost function values for each noise level at the end of each simulation, and at the end of the shortest inversion, are compared.  The shortest inversion was for the highest noise level (10%) and had 13 iterations. Comparison of both plots in Figure 11 gives an idea of the quality improvement between the 13th iteration and the last iteration for each noise level. Thus, for higher noise levels, the inversion quality shows little improvement with an increasing number of iterations. The anomaly identification (inverse analysis) results for the noise levels of 1%, 5% and 10% for setup 3 are depicted in Figure 12. A visual inspection of the FWI results in Figure 12 shows that for each increasing level of noise in the measurement data, there is a corresponding deterioration in the inversion quality, especially in the dam foundation. This behavior is also visible when looking at the cut profiles for the model parameter V s in Figure 13, because the anomaly identification in the foundation part becomes worse as the noise level increases. However, due to the sensor/receiver arrangement on the dam slope, informative wavefields transmitted through the structure enable better identification of anomalies (in the V s model) in the dam structure. Considering the results obtained in Figures 11-15, it can be deduced that, depending on the noise level in the measured data, an increase in the number of iterations does not necessarily yield better results in relation to the required additional computational effort. The required computation effort (in terms of time) for 13 FWI iterations considering different noise levels is presented in Table 2, which also lists the required iterations for each of the frequency stages depending on the noise level. The analysis is computed on the VEGAS computer cluster at the Digital Bauhaus Lab. A total of 16 processors are used, with 8 processing elements along the x-direction of the finite-difference grid and 2 processing elements along the z-direction of the numerical domain.
In practice, noise interference during data acquisition cannot be totally eliminated. Thus, it becomes necessary to obtain an optimum number of iterations, which leads to acceptable inversion results, prevents 'over-fitting' and, at the same time, reduces the computational costs. The following estimate for the reconstruction error in Equation (8) motivates the use of a discrepancy principle that establishes a trade-off between the two terms on the right-hand side of where m denotes the model parameter, () * represents the true model, () n is the iteration number, δ is the noise level being considered and () δ denote quantities effected by noise. Therefore, an attempt is made to efficiently calculate the optimal number of iterations that is required to identify the anomalies in the structure, especially for noisy measurement data. This is realized by stopping it for the first time F(m) − u δ ≤ δ, where δ ≥ 1 and δ is a measure of the noise level in the (synthetically) collected data. Therefore, the reconstruction errors of the model parameters V p and V s are plotted over the number of iterations in Figure 14. Here, it can be observed that the model V p has its minimum at iteration number k = 10, except for the noise level of 10 % that stops earlier. In contrast, the minimal reconstruction error of the model V s lies at k = 6. Although the V s model provided better inversion results than V p , a comparison between the FWI results after k = 6 and 28 iterations, which is the final one, is shown in Figure 15. The comparison also intends to check whether the models with the least reconstruction error provide the best results. (b) Figure 14. Reconstruction error considering 0%, 1%, 2%, 5% and 10% noise in the data.  The results achieved in Figure 15 differ significantly from what was expected looking at Figure 14. Although the 28th iteration has a higher reconstruction error, it gives the best results according to the visual assessment (i.e., identification of the most damaged regions) at the expense of increased artifacts. A closer look at the FWI results for higher iteration numbers (k = 28) shows that, though the inversion quality is better, there exist more artifacts, indicating damages in regions that are undamaged. Thus, these artifacts result in a higher reconstruction error, as seen in Figure 14.

Cyclic Multi-Frequency Stage Inversion
To reduce the effect of these artifacts at higher iterations and stabilize the inversion, cyclic multi-frequency analysis with time damping is applied. Here, for each time frame (t lim [s]), the inversion is carried out in different stages, starting from lower to higher frequencies (i.e., from 0.4 kHz to 1 kHz) to alleviate the non-linearity and the ill-posedness of the problem and returning again to lower frequencies. Amplitudes of the acquired signal lying outside of t lim are damped to 0. The output of each time frame is used as an initial model for the next, starting the inversion again, from low to high frequency and back again, as sketched in Figure 2. One cycle is generally completed at the end of each time frame (i.e., when the inversion reverts to a lower frequency after completing a higher-frequency stage). In order to ensure a good resolution of small anomalies, it is necessary to finish with a high frequency. Therefore, it is recommended to increase the frequencies in the last stages towards the maximum used frequency, as also shown in Figure 2. Information on the frequency bands and time frames in each cycle is obtained in Table 3, i.e., the cyclic multi-frequency approach is implemented with ω = 200 Hz, and three cycles are realized. A low-pass frequency filter is applied. The corner frequencies and the window considered are indicated in Table 3. Each time window has an origin of 0 s. In opposition to the results gained without the cyclic multi-frequency stages, using these frequency cycles leads to similar identifications of the model parameters which do not strongly depend on the noise level. The results of the cyclic multi-frequency inversion approach are given in Figure 16 under the consideration of different noise levels. When comparing Figures 12 and 16, the identification of the two areal anomalies above the tunnel of the dam has become much clearer regarding all three model parameters, but they are especially visible for the parameter V s . One drawback of the cyclic multi-frequency stage analysis is the less accurate identification of the crack at the lower dam area. Nevertheless, the crack is visible, yet not as precise as in Figure 14. In total, using cyclic stages for multiple frequencies is reducing artifacts, which is particularly noticeable in the foundation. A second effect is that both model parameters, the density and the primary wave velocity, became more informative than for the multi-frequency inversion approach, as the anomalies can also be identified when looking at these parameters, even though the identification by V s is still recommended in comparison.
In addition to the above, the stability of this regularized method compared to the multifrequency stage inversion can be seen in the reconstruction error plots in Figure 17. For the V s plot (Figure 17b), it can be observed that, as the noise increases, the reconstruction plot is decreasing step-wise at lower iteration numbers, but also, it is stabilizing with higher iteration numbers. The total number of iterations is more or less the same for the noise levels of 1%, 2% and 5%, but it is increasing for 10% noise and also for the case of no noise (0%). In total, there are more iterations required for the cyclic multi-frequency approach than for the multi-frequency stage inversion, as there are more frequency bands used during the inversion process, but nevertheless, there are fewer artifacts present in the material parameter results and the quality of reconstruction is increased. This trend is barely noticeable in the V p plot, but there is some enhancement in the reconstruction error of V p after the FWI compared to the initial model. This is because most of the damages/anomalies are identified best by the S-waves (V s ). 20 40   In addition to the noise in the data, an important role is played by the starting model in identifying the material parameters correctly. This circumstance is displayed in Figure 18, where the FWI results for −10%, +2% and +10% differences in the material parameters of the starting model from the ones of the dam's as-built parameters (see Figure 4b) are shown. When the error in the starting model is rather small (e.g., +2%) the results from the FWI are comparable to the ones when there is no deviation in the starting model. Comparing the usage of −10% (Figure 18a) and +10% (Figure 18c) difference of the material parameters in the initial model, it becomes obvious that the reconstruction of the anomalies performs better in terms of underestimating the real material properties. In that case, all but the middle anomaly can still be detected using the cyclic multi-frequency approach. As the reduction in the material properties of the second anomaly is exactly 10%, it is incorporated by the whole dam. Therefore, it can be stated that anomalies which deviate less or equally from the real model than the starting model cannot be determined by the FWI. Overall, the crack in the lower part of the dam is detected for all displayed differences in the starting model. Small changes in the starting model do not worsen the identification of the real material parameters. Additionally, it is recommended to stay on the safe side with the initial model when underestimating the material properties. From a practical point of view, uncertainties about the true condition in the undamaged state of more than 10% seem unrealistic.
(a) (b) (c) Figure 18. Cyclic multi-frequency inversion anomaly identification considering certain difference in dam as-built material properties in the starting model having additionally 2% noise in the data. (a) −10% difference in starting model (b) +2% difference in starting model (c) +10% difference in starting model.

Conclusions
Our findings in this research paper show the applicability of Full Waveform Inversion (FWI) in the non-destructive testing (NDT) of superstructures as a viable option for the current difficulties and inefficiencies faced in dam monitoring. To overcome the shortcomings of the standard FWI, we apply our cyclic multi-stage approach as discussed in Section 4.3. The robustness of our approach is shown by the high-resolution reconstruction of the damaged regions, considering the varying intensity of noise used to corrupt the data. Further findings from our research to be considered are as follows: 1.
The proposed FWI formulation is capable of effectively identifying and quantifying regions of weaknesses (i.e., heterogeneity) in both the dam structure and its foundation.

2.
The dam's as-built material properties are used as a starting model for the inversion. This information is, in most cases, readily available, or can be easily estimated. If the material properties are completely unknown, it is recommended to underestimate and not to overestimate them.

3.
The efficiency of this method is influenced by the data acquisition geometry. Thus, we propose an acquisition setup which encloses regions in which critical damage is expected. 4.
The damaged regions are generally of lower velocities and smaller scale; thus, the V s model, with its shorter wavelength, resolves the anomalies better than the V p and ρ models.

5.
The multi-cyclic-stage inversion prevents the reconstruction error of the V s parameter from increasing with a higher number of iterations (see Figure 4b vs. Figure 17b).

6.
A superimposition of the V p and V s models to identify damaged regions saturated with water increases the robustness of the method by leveraging the advantages of both models. 7.
The identified damages in the dam body had a better quality than the ones in the dam foundation for data corrupted with high noise levels. Thus, we propose, where possible, an acquisition setup which favors the recording of transmitted waves.
An outlook for future work would consider excitation sources placed within the dam structure, especially with regard to the reconstruction of better V p and ρ models.