Improving crosshole ground-penetrating radar full-waveform inversion results by using progressively expanded bandwidths of the data

In the last decade, time-domain crosshole ground-penetrating radar full-waveform inversion has been applied to several different test sites and has improved the resolution and reconstruction of subsurface properties. The full-waveform inversion requires several diligent executed pre-processing steps to guarantee a successful inversion and to minimize the risk of being trapped in a local minimum. Thereby, one important aspect is the starting models of the full-waveform inversion. Generally, adequate starting models need to fulfil the half-wavelength criterion, which means that the modelled data based on the starting models need to be within half of the wavelength of the measured data in the entire investigation area. Ray-based approaches can provide such starting models, but in the presence of high contrast layers, such results do not always fulfil this criterion and need to be improved and updated. Therefore, precise and detailed data processing and a good understanding of experimental ground-penetrating radar data are necessary to avoid erroneous full-waveform inversion results. Here, we introduce a new approach, which improves the starting model problem and is able to enhance the reconstruction of the subsurface medium properties. The new approach tames the non-linearity issue caused by high contrast complex media, by applying bandpass filters with different passband ranges during the inversion to the modelled and measured ground-penetrating radar data. Thereby, these bandpass filters are considered for a certain number of iterations and are progressively expanded to the selected maximum frequency bandwidth. The resulting permittivity full-waveform inversion model is applied to update the effective source wavelet and is used as an updated starting model in the full-waveform inversion with the full bandwidth data. This full-waveform inversion is able to enhance the reconstruction of the permittivity and electrical conductivity results in contrast to the standard full-waveform inversion results. The new approach has been applied and tested on two synthetic case studies and an experimental data set. The field data were additionally compared with cone penetration test data for validation.


I N T RO D U C T I O N
Detailed and high-resolution characterization of variable saturated soil-aquifer systems is critical and full of challenges, but highly important to improve the understandings of flow and transport processes of subsurface water.Small-scale soil heterogeneities can have a significant influence on the plant water uptake process, the transport process of pollutants or greenhouse gas emissions from peatland and permafrost soils (e.g., Simmer et al., 2015;Vereecken et al., 2016;Klotzsche et al., 2019a).In addition, the precision of large-scale atmospheric models is influenced by recharging and transport of subsurface water (Becker, 2006).To enhance our understandings of such processes, high-resolution images and accurately described subsurface properties, especially for small-scale heterogeneities of aquifers, are highly important (e.g., Klotzsche et al., 2019b).Thereby, some key parameters of geophysical and hydrological investigations, such as relative dielectric permittivity, electrical conductivity, hydraulic conductivity and porosity, play an important role (e.g., Binley et al., 2015).Traditional approaches to characterize aquifer properties either capture a small spatial sampling with a high vertical but poor lateral resolution, such as drillings and core sampling (Tillmann et al., 2008), or, have an average response over a large volume with a lack of detailed characterization at smaller scale, such as pumping test and remote sensing methods (Landon et al., 2001).In the last decades, geophysical methods have been developed and applied to obtain images of the near subsurface and to improve the characterization of hydrogeological properties (Binley et al., 2015), such as seismics (e.g., Doetsch et al., 2010), electrical resistivity tomography (ERT; e.g., Coscia et al., 2011) and ground-penetrating radar (GPR; e.g., Klotzsche et al., 2018).Especially, crosshole GPR has shown a great potential to characterize aquifers and has matured because of the opportunity to derive the highest possible resolution by using high-frequency electromagnetic (EM) pulses (e.g., Paz et al., 2017), and the opportunity to simultaneously derive electromagnetic wave velocities and amplitudes.The estimated EM velocities and attenuations can be transformed into relative dielectric permittivity ε r and electrical conductivity σ (e.g., Ernst et al., 2007b;Dafflon et al., 2011), which are related for example to water and clay content, respectively (e.g., Kowalsky et al., 2005;Cassidy, 2007;Klotzsche et al., 2014).For crosshole applications, the EM pulses are emitted from a dipole-type antenna in a borehole and received by an antenna in another borehole.Within the scope of hydrogeological site characterization, several studies have shown the potential of crosshole GPR by comparing GPR permittivity results with independent measurements of soil water content (e.g., Looms et al., 2008;Kuroda et al., 2009).The good consistency illustrates crosshole GPR is able to close the gap between small-scale investigations with highresolution (e.g., coring) and large-scale zones mapping (e.g., flowmeter tests).
Traditionally, ray-based approaches based on the geometrical ray theory are used for crosshole GPR data to derive tomographic images of the subsurface, which exploit firstarrival travel times and maximum first-circle amplitudes of the GPR wave (e.g., Maurer et al., 2004;Dafflon et al., 2012).For such approaches, damping and smoothing constraints are necessary to stabilize the inversion (e.g., Holliger et al., 2001;Maurer and Musil, 2004).Thereby, most applications consider only a limited angular coverage of the rays to avoid an increasing apparent velocity for increasing ray path angles (Peterson, 2001).Ray-based approaches are often not able to resolve targets smaller than the dominant wavelength and relatively smooth images are obtained, with a resolution that scales approximately with the diameter of the first Fresnel zone.In comparison to ray-based approaches, full-waveform inversion (FWI) exploits the entire information of the data (or significant parts of it), and is therefore able to provide higher resolution images within the sub-wavelength scale.Inspired by the FWI approach that was first developed and applied in the seismic community (e.g., Tarantola, 1984;Shin and Cha, 2008;Virieux and Operto, 2009), a 2D time-domain FWI for crosshole GPR data was implemented by Ernst et al. (2007a, b) and Kuroda et al. (2007), based on solving Maxwell's equations.Meles et al. (2010) improved the method of Ernst et al. (2007a) by including the vector characteristics of the EM fields, and introduced a simultaneous update of the permittivity and electrical conductivity parameters.Klotzsche et al. (2019b) provides an overview of the current developments, applications and corresponding challenges for the 2D crosshole time-domain GPR FWI.Several applications of the FWI to crosshole data sets from different environments, ranging from various aquifers to karst environments, demonstrated that the FWI was able to provide higher resolution images than the ray-based approaches, and provided images within decimetre scale.
Generally, FWI approaches can be implemented in the time and frequency domains.Both methods have pros and cons.While frequency approaches can highly minimize the calculation costs and improve the cycle skipping problem, the time-domain approaches provide the most flexible framework to apply time windowing of arbitrary geometries (Virieux and Operto, 2009).Similar to seismics, frequency-domain FWI approaches for GPR data have also been implemented.Such approaches allow, for example, to use only a few discrete frequencies of the data and allow the possibility to implement a wider range of misfit functions, which are beneficial for frequency-dependent medium properties (e.g., Lavoué et al., 2014).Ellefsen et al. (2011) inverted measured crosshole GPR data from a laboratory tank by using a 2.5D frequencydomain FWI approach.Furthermore, a 2D frequency-domain quasi-Newton approach for multioffset GPR data was implemented and tested on a synthetic model (Lavoué et al., 2014), and on GPR data acquired at carbonate rocks (Pinard et al., 2016).Although frequency-domain FWI seems to have certain benefits compared with time-domain FWI, for most experimental GPR data, low-frequency data are missing or show a low signal-to-noise ratio.Until now, almost all successful applications to experimental data have been performed using the time-domain FWI approach for GPR data (Klotzsche et al., 2019b).
The 2D time-domain crosshole GPR FWI uses a conjugate-gradient algorithm (Polak et al., 1969) to minimize the misfit function between the measured and modelled data.This minimization can be achieved by updating the model parameters ε r and σ by calculating updating directions and corresponding step lengths for both parameters according to the conjugate-gradient algorithm.It is well known that the non-linearity problem of the forward modelling is associated mainly with multiple scattering of propagating waves (Mora, 1987;Meles et al., 2012).As with most of the inversion approaches, the FWI with conjugate-gradient method is also ill posed, which means that we can have data residuals of multiple models which fit the data equally well (e.g., Backus and Gilbert, 1968;Brittan et al., 2013).In addition, the FWI results can be trapped in local minima if certain criteria are not fulfilled.One of the most important criteria is that the forward modelled data of the starting model are within half of a wavelength of the measured data in the entire inversion domain.Normally, ray-based starting models fulfil this criterion (Virieux and Operto, 2009).However, in the presence of specific subsurface structures such as high contrast layers related to, for example, the presence of the water table or small-scale heterogeneities, the ray-based starting models can be inaccurate and adaptations are required (e.g., Klotzsche et al., 2012 and2013).An experienced user and other tools, such as the amplitude analysis approach, are needed to ensure that the measured data are fully understood before starting the FWI and in the entire domain of interest the starting models fulfil the half-wavelength criterion (Klotzsche et al., 2014;Zhou et al., 2020).
To solve the problem of being trapped in a local minimum in the FWI, and to avoid a detailed pre-processing of the data, some solutions were proposed in microwave imaging and seis-mic inversion.For example, the frequency hopping method was used in frequency domain for microwave applications (Chew and Lin, 1995;Dubois et al., 2009) and for seismic data (Pratt et al., 1998;Maurer et al., 2009).Thereby, the inversion is first performed with a low-frequency bandwidth.The results of this low-frequency inversion are then used as the next starting models for the FWI, with a progressively increased frequency bandwidth.These steps are repeated until the maximum boundary of the frequency bandwidth reaches the centre frequency of the data.One problem of this approach is that, for many field experimental data sets, not enough lowfrequency data are present or that the low-frequency data are highly contaminated with noise.Therefore, several researchers tried to generate artificial low-frequency information for seismic data, for example, by using mathematical transformations based on the angle difference identity for cosine (e.g., Wang, et al., 2019) or using the modulation signal approach introduced by Wu et al. (2014).To tame the non-linearity issue of GPR inversion, also for the time-domain FWI, Meles et al. (2012) presented an inversion scheme that incorporated a frequency-dependent effective source wavelet.Thereby, the modelled data are progressively bandwidth expanded (PBED) in an iterative process based on a stepwise increased frequency bandwidth of the source wavelet, while the observed data kept the full bandwidth.This approach worked very well for synthetic data, but for experimental data, the choice of the frequency bandwidths and the inversion parameters such as the perturbation factors hindered a successful application so far.
In this paper, we introduce a novel approach that continues the work of Meles et al. (2012), by improving the permittivity starting model and the effective source wavelet, based on progressively expanded bandwidths of the modelled and observed data (PEBDD).In contrast to the approach of Meles et al. (2012), we applied tapered bandpass filters to the effective source wavelet (used for the modelled data) and the observed data.The new FWI scheme can be divided into two steps: (1) construction of the new permittivity starting model by using PEBDD and (2) performing the FWI with the new starting model and an updated effective source wavelet using the full bandwidth data.To evaluate our new inversion scheme, we performed two synthetic case studies.One study is conducted using standard ray-based starting models and the second considers an enforced smaller starting ε r model, which provides modelled data more than half a wavelength away from the measured data.After the validation of this new approach, we used the crosshole GPR FWI with PEBDD scheme to characterize the Krauthausen test site, using four inline crosshole GPR sections.

Standard full-waveform inversion scheme
In this section, we only discuss the most important data processing and inversion steps of the standard full-waveform inversion (FWI) approach.For more detail, we refer the reader to Klotzsche et al. (2019b).One critical aspect to guarantee reliable and stable crosshole ground-penetrating radar (GPR) FWI results is to define adequate starting models.These starting models should return modelled data that are within half a wavelength (λ/2) of the measured traces in the entire inversion domain to avoid cycle clipping and trapping of the inversion in a local minimum of the misfit function (Meles et al., 2012).In addition, an effective source wavelet needs to be estimated, because the source wavelet emitted into the earth system is unknown for experimental GPR data and accounts additionally for coupling effects.Such an effective source wavelet is estimated with two steps.First, an initial source wavelet is estimated using only horizontally travelling rays.The initial wavelet only considers the shape of the wavelet and the amplitude is normalized to one.This initial wavelet is in a second step updated using the deconvolution approach (e.g., Ernst et al., 2007b;Klotzsche et al., 2010) in the frequency domain, which can be described by: and where s k represents the initial source wavelet.E syn is the forward modelled data based on the ray-based starting models and s k .G is Green's function and E obs is the observed GPR data.η D and η I are prewhitening factors that are applied to stabilize the solution and avoid dividing by zero.ˆindicates frequency domain.s k+1 describes the updated effective source wavelet with an optimized phase and amplitude.If needed, this wavelet can be updated after a certain number of FWI iterations.Note that for most experimental data applications, the effective source wavelet is based on the ray-based models.
The forward modelling is based on a 2D finite-difference time-domain (FDTD) solution that solves Maxwell's equations in the Cartesian coordinates (Ernst et al., 2007a;Meles et al., 2010).In the process of the inversion, the misfit function C (ε, σ ) is defined by the differences between the observed data E obs and model-predicted data E syn (ε, σ ): where s, r and τ are transmitters, receivers and the observation time, respectively.T denotes the transpose operator.Each of the fields is locally defined at any point of space x and time t.The multiplication with the Dirac delta δ function selects from the entire wavefield the used receiver locations and observation times.Additionally, the gradients of the misfit function with respect to permittivity ∇C ε and conductivity ∇C σ are calculated by a zero-lag cross-correlation of the synthetic wavefield with the back-propagated residual wavefield: where R S can be interpreted as the back-propagated residual wavefield in the same medium as the incident wavefield E syn .The spatial delta function δ (x − x ) in equation ( 4) corresponds to the spatial components of the gradients and reduces the inner product to a zero-lag cross-correlation in time (Meles et al., 2010).Note that two different step-lengths are considered in the process of estimating how far the permittivity and conductivity models need to be updated in the gradient direction.The medium parameters are updated to reduce the misfit function equation ( 3) by updating all gradient values in equations ( 4) and (5).In this paper, we employed the gradient normalization approach as proposed by van der Kruk et al. (2015) to minimize inversion artefacts close to the transmitter and receiver positions in the vicinity of the boreholes.
As stopping criterion for the inversion, we consider that the root-mean-squared (RMS) error between observed and modelled data changes less than 0.5% between two subsequent iterations.Besides the stopping criterion, reliable FWI results also include the absence of a remaining gradient for the final models and a match between the measured and modelled data with a correlation coefficient greater than 0.8 (e.g., Klotzsche et al., 2019b).For synthetic case studies, the mean absolute error (MAE) in the 2D domain can be described as (Zhou et al., 2021): where RM i, j and FWI i, j represent the input and the FWI model value located at the cell of i, j. k and m show the number of cells in the 2D domain along horizontal and vertical directions, respectively.

Novel progressively expanded bandwidths of the modelled and observed data full-waveform inversion scheme
Different from the approach of Meles et al. (2012) that is using the expanded bandwidths for the effective source wavelet only, we introduce a new scheme that applies the progressive increase bandwidths on both measured data and effective source wavelet.Our new updated approach proposes to tame the non-linearity issue of the time-domain FWI, which can be considered as an extension of the standard FWI procedure to improve the characterization of small-scale subsurface structures.Therefore, we apply the idea of a frequencydomain approach that uses longer wavelengths with lower frequencies in the beginning of the inversion to avoid the cycle skipping problem and that defines the bandpass filters according to the centre frequency of an effective source wavelet.First, we construct a series of bandpass filters with different high cut frequencies, while keeping the low cut frequency constant.Figure 1(a) shows an example for an effective source wavelet with a centre frequency of 65 MHz and a bandwidth of 12-140 MHz.For such a wavelet, we would select a low cut frequency of 12 MHz, while the maximum cut frequency is considered larger than the centre frequency of the effective source wavelet, which in our case is 68 MHz.To smooth the bandpass, we assigned two tapers with lengths of 12 MHz and 10 MHz for the starting and ending frequency points, respectively (Fig. 1a).These tapered bandpass filters are applied in the next step to the observed GPR data and the effective source wavelet, which results after applying the FDTD in sub-modelled data with a similar frequency spectrum as the sub-observed data (Fig. 1b, green box).The comparison of the effective source wavelets with the expanded bandwidths with the different bandpass filters shows a clear effect on the phase and amplitude of the wavelets (Fig. 2).
Second, the FWI with the progressively expanded bandwidths of the modelled data and observed data (PEBDD) is performed (Fig. 1b, green box loops 1 and 2).The first start-ing models for the FWI are based on the ray-based results Model Ray ('Ray ( ε r ) and Homo (σ )' in Fig. 1b, green box).Traditionally, for the conductivity a homogeneous starting model is used based on the results of the first cycle amplitude inversion of the data (Klotzsche et al., 2019b).Using the starting models with the sub-source wavelet and sub-observed data with the same smallest bandwidth, we perform a certain number of iterations of the FWI.Note that we chose for our studies five iterations for each bandwidth expansion n (Fig. 1b, green box loop 1).Several synthetic tests indicated that this choice is the best compromise between computational efficiency and stability of the inversion.The perturbation factors for the inversion, which are necessary to define the steplengths for the gradient approach, are kept the same as for the standard FWI.The ε r and σ FWI results after this fifth iteration are considered as the next new starting models for the next sub-data with progressively expanded bandwidth, while the maximum cut frequency is increased by 4 MHz in our case (Fig. 1b, green box loop 2).These steps are repeated until the selected maximum cut frequency is reached (Fig. 1b green box loop 1: n = n max ).Until this point all data (including source wavelet and observed data) used for the inversion are progressive bandwidth expanded.
Third, the FWI with the full bandwidth data (FBD) is calculated (Fig. 1b, red box).From these results, only the PEBDD permittivity FWI results with the maximum cut frequency are considered as new permittivity starting model in the next step (see, 'low-freq FWI ( ε r )' in Fig. 1b).Together with the conductivity starting model, which is the same as for the traditional FWI (Fig. 1b, 'Homo (σ )'), we construct the new starting models Model New .Note that tests indicated that using the conductivity results with the maximum cut frequency as starting models does not improve the final results.The new starting models Model New and the traditional effective source wavelet SW Ray are used in the next step to generate an updated effective source wavelet SW Low with the full bandwidth data using the deconvolution approach (equation (2), Fig. 1b, red box).After obtaining the updated source wavelet SW Low , an updated FWI is performed with the new starting models Model New and the wavelet SW Low .Generally, a second-updated source wavelet SW New is necessary to further improve the FWI results.Thereby, the effective source wavelet SW Low is updated to SW New using the deconvolution method with the wavelet SW Low and the permittivity model from the updated FWI results with SW Low .The final updated FWI results with FBD are solved based on the second-updated source wavelet SW New and the new starting models Model New (details can be found in Fig. 1b).Meles et al., 2012).In the first part (green box) the progressively expanded effective source wavelet and the observed data are used, in which i, n and k represent all iterations using all sub-data, the numbers of bandwidth expansions, and individual iterations of each group sub-data, respectively.n max is depending on the centre frequency of the effective source wavelet.In this case, we select n max = 13.In the second part (red box), the full bandwidth data (FBD) FWI is performed, where two effective source wavelet corrections are performed and used during FWI.The final result is indicated by the second FWI with the SW New as an effective source wavelet.

Synthetic study I: Stochastic input models and ray-based starting models
To verify the aforementioned new full-waveform inversion (FWI) scheme for improving the ground-penetrating radar (GPR) FWI results, we construct realistic synthetic models of relative dielectric permittivity and electrical conductivity using a stochastic simulation (sequential Gaussian simulation).For the simulation, existing data set of the Krauthausen aquifer test site in Germany is used to generate a facies model with certain hydrological and geophysical parameters (based on Tillmann et al 2008, and more details can be found in Zhou et al 2021).The construed aquifer consisted of a three-layered  structure similar to the Krauthausen aquifer: sandy layer between 1.0 m and 4.0 m, sandy gravel between 4.0 m and 5.5 m and below coarse gravel (left side in Fig. 4 and Fig. 8).For the unsaturated zone above water table at 2.0 m, we chose a homogeneous layer with a relative permittivity of ε r = 4.4 (not shown, same for all following inversions).To generate the realistic synthetic GPR data (called observed data), we used a source wavelet (SW Real ; dashed curves in Fig. 3a and b) based on available experimental data from the cross-section B38-31 of the Krauthausen test site (Gueting et al., 2015).Similar to the acquisition of experimental data, a semi-reciprocal acquisition set-up was used with transmitter and receiver spacing of 0.5 m and 0.1 m, respectively.The realistic synthetic crosshole GPR data are modelled with the same time-domain 2D FDTD approach as used for the FWI.Similar to experimental data applications, we defined the ε r starting model for the FWI by picking the first-arrival times of the synthetic data and performed a ray-based inversion (first column of Fig. 4a; e.g., Maurer and Musil, 2004).Similar to previous studies, we choose for a homogeneous σ starting model of 13 mS/m (not shown) that was defined by testing various different homoge-neous models and used the results of the first cycle amplitude inversion of the experimental data.
To keep it as realistic as for experimental data applications, we used the ray-based starting models to estimate an effective source wavelet (SW Ray ; blue curves in Fig. 3a and b) and performed the standard FWI (first column of Fig. 4b and c).As expected, the FWI results show higher resolution images as the ray-based results within decimetre-scale resolution.Because of the known input models, we can calculate the mean absolute error (MAE) according to equation ( 6) between the resolved FWI models and the true input models based on the stochastic simulation (shown in the first column titles of Fig. 4d and e).The final results were obtained after 28 iterations, and the inverted results fulfilled the criteria for a reliable inversion.The corresponding RMS curve that is calculated between the observed and modelled radar data is indicated by the blue curve in Fig. 5 with the final RMS of 6.9945 × 10 −7 .Generally, most of the distinct features of the stochastic input models are resolved, and a good fit between the observed and the modelled data is obtained (see, Fig. 6a and b).By comparing the differences between observed data and modelled data for one exemplary data set (Fig. 6c), we find most of the regions in the tomograms a small misfit is visible, but for some domains a increased difference can still be noticed.Considering that the real effective source wavelet has a centre frequency of 69 MHz and the stochastic model has an approximate average permittivity value of 18, the corresponding wavelength of the GPR signal is 1.03 m.Comparing the model cell size of 0.09 m with the wavelength 1.03 m, it is hard to match all the features of the stochastic model especially when the contrast is relatively high.
Following our new introduced approach of the progressively expanded bandwidths of the modelled and observed data (PEBDD), we applied the different bandpass filters as shown in Fig. 1(a) to the standard effective source wavelet SW Ray (different effective source wavelets shown in Fig. 2) and the observed data.For this study, we choose the optimal five iterations for the sub-data FWI by comparing with other iterations values.The first sub-data FWI started from the filtered sub-source wavelet and filtered sub-observed data with bandwidth 12-16 MHz and the ray-based starting models Model Ray .Every five iterations the bandpass is increased by 4 MHz until the final cut-off frequency of 68 MHz is reached.The final permittivity FWI results using the maximum bandwidth sub-data are shown in Fig. 4a (title is 'Updated -ε r ').Combining the conductivity starting model with a homogeneous value of 13 mS/m, we construct the new start-ing models Model New for the following FWI results and the updated source wavelet.Following the flowchart (Fig. 1b, red box), the effective source wavelet SW Low is updated using the new starting models Model New and the wavelet SW Ray .The effective source wavelet SW Low together with the starting models Model New is used to calculate the FWI as shown in Fig. 4 (a and b) (second column) using the full data set.By comparing the RMS curve behaviour of this inversion (black graph in Fig. 5; black graph is covered by red graph before 71 iterations), we can notice that first the RMS is stepwise increased until the full data are used and the RMS curve with the full data are decreased after 23 iterations to a final value of 6.9749 × 10 −7 .The misfits between the input and resolved tomograms (Fig. 4d and e) and the data misfit (Fig. 6c) are slightly better than using the standard FWI approach.
These updated FWI results are in the last step considered to update the effective source wavelet a last time (SW New in Fig. 3a and b, red source wavelet) by using the new FWI ε r results with SW Low in the deconvolution approach.The corresponding final FWI results (Fig. 4, third column) are derived with the SW New and Model New as the starting models.Note that we also tested the new approach with the both permittivity and conductivity starting models based on the maximum bandwidth sub-data (  and Model New as the starting models; therefore, we consider in further applications only the permittivity updated starting model.For comparisons, we performed the approach of Meles et al. (2012) that only uses the sub-source wavelets in the FWI, which means using the progressively bandwidth expanded modelled data (PBED), while the observed data include the full bandwidth (Fig. 4, fourth column).Note that the generated FWI results (including ε r and σ ) with the maximum bandwidth sub-source wavelet are the new starting models (only show the ε r starting model with the title M − ε r in Fig. 4a) for the following FWI with full bandwidth data.The FWI results using the SW New show the smallest misfit of the resolved tomograms and the smallest final RMS, indicating that this approach is able to resolve the input tomograms best.Note that we are able to improve the FWI results for the permittivity and the conductivity by 13.7% and 7.7% in comparison to the standard FWI models using SW Ray by comparing the MAE values of the two approaches, respectively (Table 1).Especially, the small-scale structures close to the boreholes are clearer and more accurately resolved as by the standard method.By comparing the final RMS curves (only show the max iterations until 101 or 31; the same for Fig. 9) of the five different FWI results (green graph based on results shown in the Appendix), the second-updated FWI results with the wavelet SW New indicate the best convergence behaviour resulting in the smallest residual value of 2.8996 × 10 −7 after 50 iterations (red curve in Fig. 5).Comparisons of four different FWI ε r and σ results in 2D domain (Fig. 4b  and c) show the FWI results with the approach of Meles et al. (2012) are unrealistic.By computing the absolute errors between these different FWI results and the real input models in 2D domain (Fig. 4d and e), we can find that the secondupdated FWI results with SW New are close to the input models indicated by the smallest MAE values (Table 1).By investigating the fit between the observed and the FWI modelled GPR data, generally a good fit can be observed for all the FWI results (one example shown for transmitter depth 5.69 m in Fig. 6a and b), but analysing the differences (Fig. 6c) in more detail, we notice that the FWI using SW New shows the smallest misfit.
In order to describe the regional differences of different FWI models, we computed the MAE between the input stochastic models and the FWI model results along the  horizontal direction (Fig. 6d) and the vertical direction (Fig. 6e) for ε r and σ , respectively.The smallest MAE of the horizontal direction can be observed in the central part of the tomograms between 3.0 m and 5.0 m (Fig. 6d).For all results, the MAE of the horizontal direction increases towards the boundaries of the inversions domain, meaning the boreholes of the cross-section (Fig. 6d).The MAE in the horizontal direction and an increase in MAE towards the boreholes are related to the acquisition geometry and the related distribution of the ray coverage between the boreholes.Oberröhrmann et al. (2013) showed already that the resolution is highly affected by the acquisition geometry and depends on the ray coverage.The MAE along the vertical direction shows more fluctuations around 2, while the results in horizontal direction are smoother.Interestingly, all the FWI results (except the approach based on Meles et al., 2012) show improved results with a smaller MAE in the vertical direction between 4.0 m and 5.5 m depth in comparison to other depths, where the small-scale structures are located.A small increase of the MAE at 5.5 m depth can be noticed at the lower boundary of the small-scale high permittivity zone.For both horizontal and vertical directions, the FWI results based on the SW New show the lowest MAE, in contrast to the other FWI results.Comparing the FWI results and the related MAE, we can conclude that our newly updated PEBDD scheme is effective to enhance the complicated synthetic models' FWI results for permittivity and conductivity, in contrast to the standard FWI approach.

Synthetic case study II: Permittivity starting model beyond the half-wavelength criteria
In the presence of high contrast layers, ray-based results are often erroneous and need to be updated to fulfill the halfwavelength starting model criteria for the standard FWI.Here, we want to demonstrate the potential of the PEBDD scheme, which also allows starting models that are beyond the halfwavelength criteria.Therefore, we perform a second synthetic case study with the same stochastic input models as before and enforce the permittivity starting model to provide modelled data more than half a wavelength away from the measured data by using a smaller ε r model.The changed ε r starting model was obtained by subtracting three from the normal ray-based ε r model in the entire domain (Fig. 8a, title is 'Ray -ε r (−3)').The conductivity starting model is unchanged with a homogeneous value of 13 mS/m.We follow the same approach as in the synthetic study case I and obtained the different updated effective source wavelets (Fig. 7) and FWI results (Fig. 8).First, we estimate the effective source wavelet SW Ray (−3) based on the erroneous starting models Model Ray−3 .We can clearly notice that the SW Ray (−3) based on the erroneous starting models is shifted in time to the right (Fig. 7a).This is indicating that the permittivity starting model is currently too far away from the input models and that the wavelet is compen-sating for this by starting later in time (blue curves in Fig. 7a  and b).Note that a good effective source wavelet needs to start at 0 ns (Klotzsche et al., 2019b).The FWI results (first column in Fig. 8b and c) based on these starting models and the effective source wavelet SW Ray (−3) fulfil the stopping criteria after 24 iterations and no remaining gradients are present.The FWI RMS curve by using the source wavelet SW Ray (−3) is shown with the blue curve in Fig. 9. Generally, the resulting FWI results show a lower ε r than the input model, demonstrating that the ε r result gets trapped in a local minimum of the inversion process (Fig. 8b, first column).This is also indicated by the differences between the input model and FWI results in Fig. 8(d), where many regions are present that show differences with more than three in permittivity.Note that this inversion still fulfils most of the criteria for a good inversion process.The only indicators that the inversion is too far away from the input model are provided by the effective source wavelet that is shifted in time and the high MAE (not available for experimental data applications).To fulfil the criteria that the effective source wavelet needs to start at 0 ns, we shifted in the next step the source wavelet by −4 ns in the time domain SW Ray (−3shift) (cyan curves in Fig. 7a and b) to ensure that it starts at 0 ns.The corresponding FWI results (Fig. 8b and c, second column) show a better reconstruction of the permittivity results than the previous inversion, and the final RMS value after 30 iterations is 8.3326 × 10 −7 (cyan curve in Fig. 9).Meanwhile the FWI σ results are better than the FWI σ results with SW Ray (−3) (Fig. 8c, columns first and second).
Similar to the previous synthetic case І, we apply the same bandpass filters for the same number of iterations to the shifted effective source wavelet SW Ray (−3shift) to enhance the reconstruction of the FWI results with the PEBDD approach.For the first sub-data FWI, the first starting models Model Ray−3 are used.The final FWI ε r results using the maximum bandwidth sub-data (Fig. 8a, updated -ε r (−3)) and the homogeneous σ model with 13 mS/m are used as the new updated starting models Model New−3 .It can be noticed that the updated permittivity starting model shows generally higher values in the entire domain and the model is much closer to the ray-based model used in the synthetic study I, indicating that the PEBDD approach is able to enhance the erroneous starting model (compare Figs 4a and 8a).The updated starting models are considered to derive SW Low before performing the FWI.
The updated FWI results using SW Low show improved permittivity results and more continuous structures, although the conductivity performs less good indicated by the image differences (Fig. 8b and c, third column) and the higher RMS value for the final iteration (Fig. 9, black graph).The final permittivity results of this inversion and SW Low are used to update the effective wavelet a last time to SW New (Fig. 7a and b, red source wavelet).The corresponding FWI results (Fig. 8b and c, Figure 13 Standard FWI results of (a) ε r and (c) σ of the four cross-sections.Updated FWI results of (b) ε r and (d) σ based on the new ε r starting models (Fig. 11b) and corresponding updated new effective source wavelets (solid graphs in Fig. 12).Dashed lines mark the CPT data locations between the boreholes.fourth column) after 50 iterations are derived with the SW New and Model New−3 as the starting models.These results show an improved reconstruction of the permittivity and the conductivity results in comparison to the previous FWI results.Furthermore, comparisons of the absolute error tomograms for different FWI ε r and σ results in 2D domain (Fig. 8d and e) show that the second-updated FWI results with the wavelet SW New are the most accurate reconstruction of the input models.
The MAE values (shown in titles of Fig. 8d and e) and behaviour of the RMS curves (Fig. 9) are best for the secondupdated FWI results with the wavelet SW New and indicate the best FWI results (more details in Table 2).Analysing the misfit between the observed and FWI modelled data shows that the second-updated FWI results provide the best fit and indicates that this inversion obtained a model that describe the data well and best, while for the other three FWI results a significant misfit can be observed (Fig. 10a and b, exemplary  for transmitter at 5.69 m depth).To describe the regional differences of different FWI models, we compute the MAE between the stochastic input models and the FWI model results along the horizontal direction (Fig. 10c) and the vertical direction (Fig. 10d) for ε r and σ , respectively.Thereby, we can notice that the FWI ε r results with SW Ray (−3) display the largest MAE values for both directions.The FWI σ results with SW Low (black curves) have a large MAE, which implies that they cannot be resolved by only using the first-updated effective source wavelet in time when the starting models are too far away from the real models.Finally, we can conclude that the progressively expanded bandwidth scheme can not only improve the FWI results, but also that it is able to retrieve accurate FWI results for starting models more than a half-wavelength away from the measured data.Therefore, a lot of previous detailed work to construct good starting mod-els for experimental GPR data can be reduced and the application to field data could be much easier.

Experimental ground-penetrating radar data studies
As a final test, we applied the new PEBDD approach to an experimental data set from the Krauthausen test site in Germany.The Krauthausen study site is located approximately 10 km northwest of the city Düren in Germany and a detailed description of the site can be found in Vereecken et al. (2000).In the last decades, many hydrological and geophysical field techniques have been applied on this site to study the aquifer spatial distribution and flow characteristics, including flowmeter tests (Li et al., 2008), tracer experiments (Vereecken et al., 2000;Vanderborght and Vereecken, 2001), cone penetration tests (Tillmann et al., 2008) and GPR measurements (e.g., Gueting et al., 2015).Here, we utilized the measured crosshole GPR data of Gueting et al. (2015) and Zhou et al. (2021) using 200 MHz antennae (Sensors and Software) to test the PEBDD FWI scheme for four cross-sections between five boreholes.The GPR data were acquired using a semi-reciprocal acquisition set-up with a transmitter and receiver spacing of 0.5 m and 0.1 m, respectively.Before the FWI of experimental GPR data can be performed, some pre-processing steps are necessary (more details can be found in Klotzsche et al., 2019b).After applying a standard dewow filter and defining the borehole coordinates for the antennae positions, we need to determine an accurate time zero of the received radar signals.Here, we employed an improved time-zero correction method based on the ZOP-MOG cross-correlation method (Oberröhrmann et al., 2013).The ray-based method was applied for the four GPR crosssections to obtain the ray-based permittivity starting models (Fig. 11a).Similar to previous studies, the electric conductivity starting models were chosen to be homogeneous with 13 mS/m in the entire domain (not shown).In addition, to reduce the influence of 3D wave propagation phenomena, the 3D data need to be transformed into 2D according to the approach by Bleistein (1986).The effective source wavelet SW Ray (Fig. 12, dashed curves) is estimated based on the ray-based starting models according to the standard procedure.Using the wavelet SW Ray and the corresponding starting models, the standard FWI is performed (Fig. 13a and c).Note that in contrast to Gueting et al. (2015) the time-zero correction of the GPR data was improved and corrected (more details in the Corrigendum; Gueting et al., 2020).Therefore, the final standard FWI results show generally higher permittivity of about 3-4 and lower electrical conductivity results, while the structures are similar.
Similar to the previous synthetic case studies, we filtered the effective source wavelets and corresponding measured GPR data by using the PEBDD scheme.The smallest filtered sub-data bandwidths are 12-16 MHz for four crosssections.The bandpass widths are increased every five iterations by 4 MHz until the maximum bandwidth sub-data are reached.Note that we selected different highest cut corner frequencies according to centre frequencies of the different crosssection data sets.For the cross-sections B38-31 and B62-30, we used 68 MHz, and for the B32-38 and B31-62 data set, we applied 72 MHz.The final FWI ε r results using the maximum bandwidth sub-data (Fig. 11b) and the homogeneous σ model were considered as the updated starting models and applied in a following FWI under the full bandwidth data.Following the flowchart in Fig. 1(b), the effective source wavelet SW Low is updated using the new starting models and the effective source wavelet SW Ray .The updated source wavelets SW Low together with the new starting models are used to calculate the FWI results (not shown).Finally, the second-updated effective source wavelets SW New are derived using the new FWI ε r results and SW Low in the deconvolution approach (Fig. 12, solid curves).The corresponding updated final FWI results (Fig. 13b and d) show very similar structures as the standard FWI results, but the structures are more continuous closer to the boreholes.The RMS values between the observed and modelled data based on selected different FWI models show that the updated FWI results are smaller than those of the standard FWI (Table 3), indicating improved results of the updated FWI.
To validate the different FWI results, we compared them with porosity data based on cone-penetrations test (CPT) measurements (Yang et al., 2013;Zhou et al., 2021).Note that the CPT data were acquired mostly in the centre of the crosssections.First of all, we converted the porosity values φ of the CPT to the relatively permittivity using the three phase complex refractive index model (CRIM) in saturated zone (e.g., Gueting et al., 2015), where ε s and ε f are the permittivity of the solid and the fluid, respectively.
Considering a water temperature of 10°C we use a ε f of 84 and we apply a ε s value with 4.5 based on literature values of quartz (e.g., Eisenberg and Kauzmann et al., 2005).We compared the two FWI ε r results with the CPT data along the vertical profile locations (dashed lines in Fig. 13).Comparisons of these 1D vertical results indicate that the updated FWI ε r results (Fig. 14, red lines) are closer to the CPT data (Fig. 14, black lines) than the standard FWI ε r results (Fig. 14, blue lines).Quantitative comparisons are described in Table 3, where we computed RMS and correlation coefficient R with larger than 0.8 between the CPT and two FWI results.In general, the updated FWI results with the PEBDD scheme are better than the standard FWI results, indicating that the new PEBDD approach improved the experimental GPR FWI results.Note that the RMS value (for 1D CPT) of the updated FWI for boreholes 38-31 is larger than the standard FWI (green value in Table 3).A possible explanation could be that the standard FWI was already very good because many studies have been already performed with this data set, and many optimizations have been performed with the standard approach.For the other data sets, less intense tests have been applied for the standard FWI.

C O N C L U S I O N S A N D O U T L O O K
We have presented a new crosshole ground-penetrating radar (GPR) full-waveform inversion (FWI) approach using progressively expanded bandwidths of the measured and modelled data (PEBDD) with frequency bandpass filters to improve the reconstruction of the subsurface properties.The new PEBDD scheme was applied to two synthetic case studies and an experimental data set from the Krauthausen test site in Germany.Thereby, we have demonstrated that the new scheme is able to improve the FWI results by taming the problem of the inversion becoming easily trapped in a local minima caused by non-linear and ill-posed problems.The introduced PEBDD scheme applies designed bandpass filters to the effective source wavelet and the observed data.Thereby, the different sub-data sets with various bandwidths are progressively expanded after the centre frequency of the data is reached.The FWI results of this PEBDD provide an updated ε r starting model, which is used to update an effective source wavelet in the deconvolution approach and to perform the FWI with full bandwidth data.Further, we update the effective source wavelet SW New .Considering this updated effective source wavelet and the starting models based on the PEBDD scheme, the resulting FWI results show a better reconstruction of the medium properties and the final root-mean-square values are decreased in contrast to the standard FWI.
Two synthetic cases have proven that the PEBDD scheme is robust and can also reconstruct the permittivity and electrical conductivity results for starting models that are more than half a wavelength away from the measured data.To fulfil the starting model criteria, the modelled data need to be within half a wavelength of the measured data and normally need an accurate processing of the data and a detailed understanding of the trained user.Problems arising from erroneous starting models are sometimes hard to notice, and therefore an improved approach such as the PEBDD can significantly improve the applicability of the FWI.For the experimental GPR data, we have computed the standard and the updated FWI results using the PEBDD scheme for four cross-sections.We compared the different FWI results with CPT data from the centre of the cross-sections, where the updated FWI results provided the best correlation to the CPT data.Finally, we concluded that not only the progressively expanded bandwidth approach can improve the FWI results, but it is also able to retrieve more accurate FWI results for starting models more than a half wavelength away from the measured data.Therefore, a lot of previous detailed work to construct good starting models for experimental GPR data can be reduced and the application to field data will be much easier.In future work, we will also investigate the possibility to apply the new PEBDD approach to other frequencies such as for 100 MHz antennae and to different test sites.

AC K N OW L E D G E M E N T S
The first author acknowledges the financial support from China Scholarship Council (Project No. 201506340136).We gratefully acknowledge Peleg Haruzi for providing help in the process of constructing the synthetic stochastic models.Especially, we are grateful to Jessica Schmäck for her constructive suggestions.Further, we gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA at Jülich Supercomputing Centre (JSC).
Open access funding enabled and organized by Projekt DEAL.

A P P E N D I X
For the synthetic case study I, we constructed a second set of starting models Model New ( ε r + σ ), which were derived from the final FWI ε r and σ results using the maximum bandwidth sub-data (  ), it can be noticed that the new updated FWI results with SW New ( ε r + σ ) show a less good fit than the second-updated FWI results with SW New (Table 1).There-fore, we abandon this approach of using both permittivity and conductivity results derived by the FWI using the maximum bandwidth sub-data and consider for all other studies the homogeneous σ starting model with 13 mS/m. ©

Figure 1
Figure 1 (a) Example for an effective source wavelet amplitude spectrum and corresponding stepwise increased bandpass filters (shown schematically by the black and red horizontal bars).These filters are progressively expanded after the centre frequency is reached and applied to the observed data and the effective source wavelet.The bandpass filters start at the lowest frequency of 12 MHz (taper length 12 MHz) indicated by a vertical dashed line.The high cut frequency is stepwise expanded every five iterations of the FWI and has a taper of 10 MHz.After the highest corner frequency of 68 MHz is reached (centre frequency of 65 MHz), all subsequent iterations use the full bandwidth of an updated source wavelet and the observed data.(b) Flow chart of the new proposed approach (adapted fromMeles et al., 2012).In the first part (green box) the progressively expanded effective source wavelet and the observed data are used, in which i, n and k represent all iterations using all sub-data, the numbers of bandwidth expansions, and individual iterations of each group sub-data, respectively.n max is depending on the centre frequency of the effective source wavelet.In this case, we select n max = 13.In the second part (red box), the full bandwidth data (FBD) FWI is performed, where two effective source wavelet corrections are performed and used during FWI.The final result is indicated by the second FWI with the SW New as an effective source wavelet.

Figure 2
Figure 2 Comparisons of the selected sub-source wavelets with the different bandpass filters of the new PBEDD approach used for synthetic case study I in time domain.(a) Shows the source wavelets with the absolute amplitudes and (b) the corresponding wavelets normalized to the corresponding minimum (or maximum) for a better comparison of the phase.

Figure 3
Figure 3 Comparisons of the (a) effective source wavelets and the (b) corresponding frequency spectra of the different steps of the new PBEDD approach used for synthetic case study I.The legend values in (b) indicate centre frequencies for different source wavelets.The source wavelet used to generate the 'observed data' is indicated with a dashed black line.Note that SW New ( ε r + σ ) is the effective source wavelet based on the results shown in the Appendix.

Figure 4 Figure 5
Figure 4 Overview of the FWI results using the different approaches and starting models.Left, the real input models based on the stochastic simulation (Zhou et al., 2021).(a) Permittivity starting models and corresponding FWI results for (b) ε r and (c) σ for the different FWI approaches.The applied effective source wavelets are named in the legends.Image plots of the absolute error between the real input models and final FWI results are shown for (d) ε r and (e) σ for the different approaches.The mean absolute errors of the entire 2D domain are shown in parentheses.
Fig. A.1a).The obtained FWI results (Fig. A.1b) indicate the results with SW New ( ε r + σ ) did not improve significantly than the FWI results with the SW New

Figure 6
Figure 6 (a) Observed data, (b) modelled data based on the FWI results and (c) differences between the observed and modelled data for one exemplary data set of transmitter location at 5.69 m depth (see input models and black arrow for the transmitter location in Fig. 4).Note the amplitudes in (a-c) are normalized to the maxima amplitude of the real observed data (shown range from −7 × 10 −1 to 7 × 10 −1 ).The mean absolute errors of the permittivity and conductivity between the real input and the different FWI results are shown along (d) the horizontal cross-section and (e) vertical direction.

Figure 7 Figure 8 Figure 9
Figure 7 (a) Comparisons of the different effective source wavelets used for synthetic case study II.The source wavelet used to generate the 'observed data' is indicated with a dashed black graph.(b) Frequency spectra comparisons from 0 to 140 MHz for the five different effective source wavelets.The legend values in (b) indicate these effective centre frequencies for different steps.

Figure 10
Figure 10 (a) Modelled data based on the FWI results and (b) differences between the observed and modelled data for the transmitter at 5.69 m depth (see input models and black arrow for the transmitter location in Fig. 8).Note the amplitudes in (a) and (b) are normalized to the maxima amplitude of the real observed data (shown range from −7 × 10 −1 to 7 × 10 −1 ).The mean absolute errors of the permittivity and conductivity between the input models and the different FWI results are shown along (c) horizontal cross-section and (d) vertical direction.

Figure 11
Figure 11 (a) Ray-based ε r results and (b) FWI ε r results using the selected maximum frequency bandwidth sub-data.Note that results shown in (b) are used as starting models for the full bandwidth inversion of the experimental Krauthausen data.And a homogeneous conductivity starting model of 13 mS/m was used for all inversions.

Figure 12
Figure 12 (a) Comparisons of the effective source wavelets estimated based on ray-based models (dashed graphs) and the updated source wavelets of the PEBBD scheme (solid graphs) for the four crosssections in time domain for the experimental data of the Krauthausen test site.(b) Corresponding frequency spectra of standard and updated effective source wavelets.The legend values in (b) indicate these effective centre frequencies for different cross-sections.

Figure 14
Figure 14 Permittivity comparisons derived from cone penetration test (CPT) data (black), standard FWI ε r results (blue) and the updated new FWI ε r results (red) along CPT profiles between boreholes (see Fig. 13 for locations).
Fig. A.1a).Following the flowchart (Fig.1b), an effective source wavelet SW Low ( ε r + σ ) (not shown) is updated using the new starting models Model New ( ε r + σ ) and the wavelet SW Ray .This updated wavelet SW Low ( ε r + σ ) to-gether with the starting models Model New ( ε r + σ ) are used in the next step to calculate the FWI results using the full data set (not shown).These updated FWI results (including ε r and σ ) are in the last step considered to update the effective source wavelet SW New ( ε r + σ ) (Fig.3a and b, green source wavelet).The corresponding final FWI results (Fig. A.1b) are derived with the SW New ( ε r + σ ) and Model New ( ε r + σ ) as the starting models.By comparing the final RMS curves after 45 iterations (green graph in Fig.5), the second-updated FWI results with the wavelet SW New (RMS = 2.8996 × 10 −7 ; red curve) indicate a better convergence behaviour and a smaller RMS than the new updated results with the SW New ( ε r + σ ) (RMS = 3.8620 × 10 −7 ; green curve).Computing the absolute

Figure A. 1
Figure A.1 (a) The permittivity and conductivity starting models based on the FWI results using the maximum expended bandwidth sub-data.(b) Corresponding FWI results using the updated starting models and the updated effective source wavelet.Effective source wavelet is named in titles.(c) Image plots of the absolute error between the stochastic input models and the final FWI ε r and σ results.The mean absolute errors MAE of the entire 2D domain are shown in parentheses of titles.

Table 1
Comparison of the different FWI approaches for synthetic case study I using the mean absolute error between the real input models and the different FWI results for the entire 2D domain, and the root-mean-squared RMS error between observed and modelled radar traces represents residual values.Percentages in parentheses indicate the ratio of the single FWI RMS to the standard FWI RMS with SW Ray , while a decreased value means the higher improvement efficiency

Table 2
Comparison of the different FWI approaches of the mean absolute error between the real input models and the different FWI models for the entire 2D domain for synthetic case study II.RMS represents residual values between observed and modelled radar traces.Percentages in parentheses indicate the ratio of the other FWI RMS to the standard FWI RMS with SW Ray (−3) , the lower value means the higher improvement efficiency

Table 3
Comparison of the different FWI approaches for the experimental data set of the Krauthausen test site.RMS represents residual values between observed and modelled radar traces.Percentages in parentheses indicate the ratio of the new FWI RMS to the standard FWI RMS with SW Ray , the lower value means the higher improvement efficiency.Correlation coefficient R and RMS between 1D ε r FWI and the CPT data (CPT porosity has been converted into ε r ) represent the reliability of the FWI results, the larger R and the lower RMS values mean the FWI results are closer to the real values