Work f l ow to numerically reproduce laboratory ultrasonic datasets

A.Biryukov,N.Tisato,G.Grasselli

Department of Civil Engineering,University of Toronto,Toronto,Ontario,Canada

Work f l ow to numerically reproduce laboratory ultrasonic datasets

A.Biryukov*,N.Tisato,G.Grasselli

Department of Civil Engineering,University of Toronto,Toronto,Ontario,Canada

A R T I C L E I N F O

Article history:

Received 10 September 2014

Received in revised form

9 October 2014

Accepted 10 October 2014

Available online 29 October 2014

Numerical methods

Calibration

Velocity model

Bentonite

Viscoelastic wave propagation

The risks and uncertainties related to the storage of high-level radioactive waste(HLRW)can be reduced thanks to focused studies and investigations.HLRWs are going to be placed in deep geological repositories,enveloped in an engineered bentonite barrier,whose physical conditions are subjected to change throughout the lifespan of the infrastructure.Seismic tomography can be employed to monitor its physical state and integrity.The design of the seismic monitoring system can be optimized via conducting and analyzing numerical simulations of wave propagation in representative repository geometry. However,the quality of the numerical results relies on their initial calibration.The main aim of this paper is to provide a work f l ow to calibrate numerical tools employing laboratory ultrasonic datasets.The f i nite difference code SOFI2D was employed to model ultrasonic waves propagating through a laboratory sample.Speci f i cally,the input velocity model was calibrated to achieve a best match between experimental and numerical ultrasonic traces.Likely due to the imperfections of the contact surfaces,the resultant velocities of P-and S-wave propagation tend to be noticeably lower than those a priori assigned.Then,the calibrated model was employed to estimate the attenuation in a montmorillonite sample.The obtained low quality factors(Q)suggest that pronounced inelastic behavior of the clay has to be taken into account in geophysical modeling and analysis.Consequently,this contribution should be considered as a f i rst step towards the creation of a numerical tool to evaluate wave propagation in nuclear waste repositories.

?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

Isolation of high-level radioactive waste(HLRW)is an important issue which must be thoroughly addressed.HLRW is generally enveloped in multiple engineered and natural barriers and placed in deep geological repositories.Such a technique is currently considered as a viable and reliable option to safely isolate HLRW from the aquifers and the biosphere(Chapman and McCombie, 2003;Alexander and McKinley,2007).Montmorillonite is a swelling clay that is extensively used as the base material for those engineered barriers as it acts as an impermeable(hydraulic conductivity k≈10-14m/s)seal between the HLRW and the host rock as it saturates(Lajudie et al.,1994).

In the vicinity of the HLRW containers,the temperature,pressure,and water content of the barrier are expected to increase dramatically over few years after completion.The reasons for that are(a)radioactive decay of HLRW,(b)swelling of the clay and(c) imbibition from the surrounding aquifer,respectively(Villar et al., 2005;Alonso et al.,2008).Therefore,the sealing material will change the physical state and its integrity can be also jeopardized. This might lead to harmful leakage of noxious condensate into host rock and surrounding aquifers.

The aforementioned physical changes affect the elastic and viscoelastic properties of the montmorillonite,such as longitudinal and shear wave velocities(VPand VS,respectively),density(ρ),and seismic attenuation.Tisato and Marelli(2013)showed that a variation in con f i ning pressure(pc)ranging between 0 MPa and 20 MPa induces the increases in VPand VSup to 90%approximately.Similarly,the increase of water saturation(wc)from 10%to 52%at pc＜10 MPa causes increases in VPand VSup to 50%.These changes in velocities will result in a variation of the seismic signal transmitted through the montmorillonite barrier.Therefore,once the variation of the elastic parameters due to the changes in physical conditions is known,seismic monitoring may be used as a nonintrusive tool to track the condition of the plug.

The equipment employed in seismic monitoring consists of a set of emitters and receivers,installed in the proximity of the bentonite barrier(Manukyan et al.,2012).The input signal is periodically sent by the emitters,propagates through the barrier,and is recorded by the receivers.First-arrival times are evaluated based on the recorded waveforms for each cycle.If pcor wcof montmorillonite changes between any two cycles,the corresponding change in theelastic parameters will be illustrated by a f i rst-arrival time shift.As soon as changes in f i rst arrivals are detectable,full waveform analysis may be applied to locate the areas of anomalous pressure or water saturation(Manukyan et al.,2012).However,high-quality and well-calibrated data acquisition systems and experimental repeatability are required to correctly re f l ect the variation in water saturation or pressure in the barrier.Therefore,the numerical simulation of the procedure in complex repository geometry should be conducted and thoroughly analyzed to aid in design and optimization of the monitoring system(Marelli et al.,2010).

Numerical tools provide accurate results if calibrated with rigorous laboratory investigation.To the best of our knowledge, only Saenger et al.(2014)incorporated the calibration of numerical tools to support laboratory experiments(and vice versa)and to simplify the interpretation of the obtained results.Often,the analytical prediction of the recorded waveform across the geometry of the experiment investigating elastic properties of a certain material is cumbersome.Saenger et al.(2014)simulated the ultrasonic wave propagation velocity measurements in the rock samples tested in a Paterson gas-medium apparatus.It is shown that dispersive and wave conversion effects caused by the presence of assemble items(e.g.jackets or buffer rods)may hinder the f i rstarrival time.As a result,the estimated values of VPand VSwill be erroneous and deviate from values a priori assigned in the numerical model.Thus,numerical modeling of the commonly used experimental techniques is of paramount importance for accurate interpretation of the results.

The purpose of this particular contribution is twofold:(a)to propose a methodology to calibrate 2D(two-dimensional)numerical tools to yield synthetic traces,accurately reproducing the experimental results of ultrasonic wave propagation in the material of interest;(b)to use the calibrated framework and employ the iterative optimization technique as an attempt to evaluate the variation of attenuation in bentonite caused by the water saturation.The results of the calibration are validated by the comparison of the synthetic seismograms of shear wave propagating in bentonite with the experimental dataset reported in Tisato and Marelli(2013).Eventually,the long-term scope of our research is to create a numerical tool to simulate the seismic monitoring of physical changes in a generic HLRW barrier.

2.Methods

In this section we explain the strategy that allows us to numerically reproduce the experimental ultrasonic traces in montmorillonite samples collected in the laboratory by Tisato and Marelli(2013).

2.1.Laboratory setup

Tisato and Marelli(2013)determined the longitudinal and transverse ultrasonic wave propagation velocities in montmorillonite,employing an axially loaded assembly comprising two aluminum caps that enclose the sample and two ultrasonic piezoceramic transducers(ultrasonic facility).The facility allowed the recording of ultrasonic signals transmitted through the sample at set temperatures while measuring the sample length(Fig.1).The setup was designed to reproduce those conditions expected at a HLRW waste repository throughout its lifespan as well as to simulate anomalous conditions outside the predicted range.

Longitudinal and transverse vibrations with a fundamental frequency of 100 kHz(Fig.1d)were generated by means of 1 MHz corner frequency compression and planar shear piezoceramic transducers,respectively.

2.2.Numerical setup

Fig.1.(a)Schematic diagram of the ultrasonic facility employed in the experiments to measure VPand VS;(b)a subdomain used for iterative calibration;(c)a photo of the ultrasonic facility in a fully assembled state;and(d)normalized signal sent to the emitter.

We performed simulations on a 2D numerical model representing the laboratory setup to compare numerical results with laboratory data.The mesh was produced following a schematicdiagram of the longitudinal cross-section of the laboratory setup (Fig.1a).The schematic section was generated automatically by a customizable script written in Matlab.The user can specify the dimensions and the sample material.Nodal properties(VP,VS,ρand quality factors QP,QS)were assigned based on the material that a node corresponded to.The viscoelastic parameters for the materials used are reported in Table 1(Boyer and Gall,1984;Auerkari, 1996;Lakes,2009).

The transfer function of the employed transducers linearly converts the voltage applied to the transducers into acceleration within the used bandwidth(1 kHz＜f＜1 MHz).We measured the voltage signal sent to the emitter using a digital oscilloscope(Fig.1c and d).Consequently,the normalized signal was applied as a body force in a horizontal or vertical direction to excite transversal or longitudinal wave,respectively.As P-and S-wave receivers and emitters have different compositions,their locations are shown in Fig.1b with a green and a red line,respectively.

The standard staggered grid(SSG)f i nite difference(FD)code SOFI2D was used to solve the viscoelastic wave propagation problem in the given geometry.SOFI2D is an open-source 2D time domain massive parallel modeling code for P-and SV-waves.There are reasons that favor our choice of using 2D numerical model over 3D(three-dimensional)model,such as(a)computational time and (b)viscoelastic model limitations.Currently,the 3D version of SOFI software(SOFI3D)is only capable of using a standard linear solid model,whereas in SOFI2D a generalized Maxwell model can be implemented.The ultimate goal of our research is to study montmorillonite,which has a pronounced viscoelastic behavior not necessarily showing a single relaxation time.Thus SOFI2D was chosen for those reasons.Speci fi cally,we employed a fourth-order fi nite difference,time-explicit scheme to evaluate the displacement fi eld in the nodes of the mesh.The time step varied in the interval 5 10-9s?Δt?8 10-9s to simultaneously maintain numerical stability throughout calibration simulations and optimize the calculation time.

The dominant frequencies of the input signal(Fig.1d)impose a constraint on the grid spacing.This constraint,known as grid dispersion,was thoroughly studied by Alford et al.(1974).Highfrequency waves with the shortest wavelengths(λmin)will seem to slow down or even stop propagating once they get sampled less than several spatial samples per wavelength.Grid dispersion may signi f i cantly distort the results leading to serious errors in interpretation of events.To mitigate the effect,Alford et al.(1974)suggestedλmin/Δh＞6 for a fourth-order FD scheme,whereΔh is the grid spacing.Therefore,for the condition to remain valid during the whole calibration procedure,the grid spacing was set to Δh=135μm,corresponding to 30 samples per shortest wavelength observed during simulations.

During the calibration procedure,the wave propagation was simulated in a subdomain(Fig.1b)of the entire experimental setting,containing the sample assembly and the transducers (“numerical ultrasonic facility”).Numerical simulation in the subdomain decreased dramatically the computational time and did not affect the f i nal results as shown in Section 3.

Table 1Material parameters used in numerical simulations.Seismic quality factors QPand QSare de f i ned following O’Connell and Budiansky(1978)as Q=Mr/Mi,where Mrand Miare the real and imaginary parts,respectively,of the complex modulus.

2.3.Reasons for calibration

The calibration procedure presented herein is an iterative process that involves running a series of wave propagation simulations in standard samples tested in the numerical ultrasonic facility (Fig.1a and b).The input signal,traveling between the emitter and the receiver,is subjected to multiple re f l ections,refractions and scattering.The simulated propagation of the input signal in the absence of a sample(“no-sample”case)was simulated to justify the need of the calibration.The discrepancies in the f i rst-arrival times, power and phase spectra between numerical and experimental ultrasonic traces for both P-and S-wave experiments were initially observed(see Section 3).The numerically calculated f i rst-arrival times tend to be lower than those observed in the laboratory. There are a number of explanations for this phenomenon that will be addressed in Section 4.Thus,the parameters of the numerical model must be calibrated to account for the imperfections of the laboratory setup and to be able to produce accurate results in both time and frequency domains.

2.4.Calibration methodology

Systematically high apparent velocities can be corrected following a combination of two solutions similar to those proposed by Madonna et al.(2012).As a f i rst option,numerical wave propagation velocities employed in the model can be effectively decreased by multiplication by a factor k [0,1],i.e.Vexp=kVnum. The second option might be to arti f i cially introduce slow wave velocity layers at the contact interfaces to account for the asperities and vary their elastic properties to achieve the desired similarity between numerical and experimental traces.The main drawback of the latter option is that the choice of the descriptive elastic parameters and the thickness is not properly constrained.Micro-CT imaging should be employed to provide the information about the geometric details and asperities roughness(Madonna et al., 2012).A decrease in arbitrarily chosen elastic parameters may lead to a spurious change in the re f l ection and transmission coef f i cients.The resultant numerical trace may have a good f i rstarrival match but is lack of mid-and long-term overlap.Therefore,we decided to implement and focus on evaluation of the effective ultrasonic velocities of the materials employed in the model.

The spatial distribution of VP,VS,andρcomprises the minimal set of medium parameters that enable the simulation of the elastic wave propagation using SOFI2D.Therefore,to in f l uence the numerical arrival times and spectral image of the waveform,one needs to vary both VPand VSof the materials independently.Thus, to solve the problem of matching numerical and experimental seismograms,it seems necessary to run a number of models,iteratively changing the elastic parameters within a certain interval.By applying a certain goodness-of-f i t(GOF)metrics to the obtained synthetics,every numerical seismogram will be assigned a certain GOF value.The minimum and maximum GOF values among the set of seismograms will indicate the best and worst f i ts with the experimental seismogram,respectively.The GOF calculation method is addressed in details in Appendix A.

For the arrivals to match,the travel time of a direct wave between the emitter and the receiver in the setup should be equal to the f i rst arrival time texpin the corresponding ultrasonic trace. Theoretically,in the case of a layered model consisting of n layers, we have

where Liand Viare the thickness and wave propagation speed in the i-th layer,respectively.As Liis f i xed by the dimensions of the setup,Eq.(1)has an in f i nite number of sets of Vi（i=1，2…，n）to satisfy it.However,the deviation of Vifrom that derived using standard methods(Vst)for the material comprising the i-th layer is expected to lie within meaningful interval,e.g.|Vi-Vst|/Vst?0.30. Moreover,the arbitrarily chosen set will change the re f l ection and transmission characteristics of the layered con f i guration.Therefore, a set of constraints should be imposed on the Viselection procedure to conserve the general wave propagation character in the laboratory framework.However,if a medium consists of N different layers,there are 2N independent variables to be varied.This fact greatly increases the computational complexity of the problem.For M different values of each parameter,there will be M2Nunique simulations.

As a result,the number of simulations needed for M unique values of KPand KSis independent of N and is M2.

One of the advantages of this relatively simple algorithm is that it allows some degree of freedom in choosing the interval of KPand KS.Therefore,the process has an iterative character:the f i rst iteration may contain only a few points spread over a large range and hence be a rough estimation.At the end of the f i rst iteration,the GOF function is calculated and its extremes are found.The successive iterations will re f i ne the vicinity of the global maximum of GOF and thus provide a better resolution for best-f i t KPand KSuntil the f i t is considered satisfying.

3.Results

3.1.Calibration of the experimental setup

Three stainless steel samples(19.84 mm,24.84 mm and 29.86 mm long)and“no-sample”models were employed to systematically analyze the discrepancy between numerical and experimental ultrasonic traces.The preliminary wave propagation simulations showed that numerical travel times of the signal were consistently lower than experimentally measured values.This fact motivated us to impose the initial range of KPand KSvalues between 0.70 and 1.05.In the f i rst set of iterations,the range was discretized into 15 equally distributed values.Consequently,225 (i.e.15 15)numerical simulations were run corresponding to all possible combinations of KPand KS.GOF values were assigned to resultant numerical traces and series of minima at KS=0.8,KP(0.8,1.0)were observed(Fig.2a).The second iteration re f i ned the vicinity of the midpoint of the series.The procedure was repeated, until the variation of GOF function was satisfactorily small across the set of traces.For illustration purposes,best-f i t results achieved in“no-sample”simulations are compared with traces corresponding to initial and underestimated velocities(Fig.2b,c and d, respectively).However,the resultant best-f i t values of KPand KSwere consistent within all four models.

There is an evident difference in the f i rst three cycles of numerical and experimental traces(Fig.2b)corresponding to slightly changed initial velocities,implying the necessity to further decrease VPand VS.This observation supports the importance of the calibration procedure to correctly reproduce ultrasonic velocity measurements.The lower boundary for the error of the calibration may be estimated as the resolution in KPand KSvalues and hence is about 3.2%(see Appendix C for details).As a result of the calibration,the initial input parameters VPand VSshould be multiplied by factors of KS=0.775?0.025,KP=0.90?0.03,respectively,for accurate modeling of the experiment.

3.2.Subdomain model upscaling

Fig.2.(a)The GOF distribution as a function of the KP,KS,and(b,c,d)numerical traces corresponding to“no-sample”model with various velocity models.Green squares and red triangles show analytically predicted P-and S-wave arrivals.

The subdomain model strategy was employed to dramatically speed up the computations.A full-scale numerical simulation with best-f i t input parameters needs to be run in order to lend validity to KPand KSevaluated via the isolated subdomain calibration. Therefore,we propagated the input signal in full-scale model with input parameters obtained through calibration and compared the resultant numerical trace with corresponding subdomainnumerical traces and experimental traces(Fig.3).The results show the satisfactory overlap between subdomain and full-size models (see Fig.3c)is well within the accuracy of the calibration procedure. Moreover,the comparison of the experimental trace and full-size synthetics exhibits an impressive overlap not only in the calibration time range(0-70μs,Fig.3a and b),but also in the following time section(75-110μs,Fig.3b).This evidence strongly supports the capability of the subdomain calibration routine to yield parameters applicable to full-size models.

3.3.Wet montmorillonite S-wave attenuation estimation

One of the aims of the current paper is to estimate ultrasonic attenuation in a wet montmorillonite sample(52%water saturation at 30.5?C).The attenuation is de f i ned following O’Connell and Budiansky(1978)as Q=Mr/Mi.The linear viscoelastic rheology employed in the simulations was based on a“generalized standard linear solid”(GSLS)model that is able to give a realistic framework and explain the observations of wave propagation through earth materials(Liu et al.,1976).

We employed a third-order GSLS model(Fig.4)to simulate a desired constant Q-spectrum within the frequency range of interest.The characteristic parameters for the spring stiffness and viscosity of the dashpot(Fig.4,inlet)were assigned using the leastsquares optimization described in details by Bohlen(2002).

The estimation of the attenuation was organized as an iterative process in a manner similar to that described above.As the velocity model of the laboratory setup has been already calibrated following the aforementioned algorithm,the parameters to be varied were only QPand QSof the sample material(i.e.montmorillonite).The range of variation was established based on the values reported in the literature.Astbury and Moore(1970)documented a quality factor QS≈5 and Karakurt(2005)obtained a similar value for QP. Therefore,we set the range of variation to Xq=[1,20]and evaluated GOF of ultrasonic traces corresponding to numerical models with QPand QSXq.Analogous to VP,VScalibration,we obtained a 2D GOF plot with a set of global minima constrained within QP[2, 13],QS[1,7](Fig.5a).The time interval involved in GOF calculation is 0-65μs.

Two viscoelastic and a quasi-elastic(Q S=300)simulation corresponding to different values of QSare shown to illustrate how variation of the quality factors affects the similarity between experimental and numerical traces(Fig.5).Best overlap is achieved at QP≈9.7,QS≈5.8 and is relatively consistent within the range of QPused.However,the quasi-elastic model shows a strong deviation from the experimental waveform and those predicted with low quality factors shortly after S-wave arrival(t=40μs,see Fig.5c). This discrepancy clearly demonstrates the limited capability of elastic numerical models to accurately reproduce ultrasonic waveform propagating through materials showing inelastic behavior.

Fig.4.An example of“constant Q-spectrum”attenuation model employed in the models corresponding to(a)Q=6 and(b)its GSLS diagram.

4.Discussion

As itis shown by Madonna et al.(2012),here we demonstrate that a velocity model requires calibration to yield accurate results.Ultrasonic velocities estimated from theory and numerical simulations using VPand VStaken from literature,are generally higher than those measured in the laboratory.As a consequence,VPand VSshould be slightly changed to mitigate the overestimation.Moreover,this variation(i.e.decrease)should preserve the original transmission and re f l ection coef f i cients,which can be achieved by varying proportionally all VPand VS.Therefore,we state that the aid of the numericalsimulations inthe interpretationofexperiments is ultimately based on the ability of the code to reproduce laboratory results.

The uncertainty in de f i ning the stiffness at the contact interfaces between the different parts of the model may be a reason for the discrepancy.FD methods are generally restricted to the treatment of continuous interfaces.In theory,the effective elastic modulus (Eeff)of a layered medium consisting of two layers(E1,E2)is

Fig.3.Numerical traces corresponding to subdomain(a)and full-size(b)models.The overlap between the traces(c)lends validity to the values of KPand KSobtained through modeling on the subdomain.

Fig.5.The GOF distribution as a function of the QPand QS(a)and numerical traces corresponding to viscoelastic(b,c)and a quasi-elastic(d,QS=300)simulations in bentonite. Green squares and red triangles show analytically predicted P-and S-wave arrivals,respectively.

In reality,layered media are in contact through asperities which decrease the effective elastic modulus of the same medium,due to a partial discontinuity,presence of the empty space,and wave scattering.For those reasons,numerical models are prone to overestimate the ultrasonic velocities.Our calibration of VPand VSwas performed using a model which represents a subdomain of the laboratory setup.Numerical simulations run on the full-size model show that this procedure remains effective until the time range used in calibration is shorter than or close to the travel time of direct waves(0-55μs in our case).If the time range of interest is longer,the re f l ection and scattering from various components of the setup has to be taken into account and is depicted in the trace (Fig.3c).The f i rst 2-3 cycles(0-45μs)of the numerical and laboratory traces are perfectly overlapped,whereas small discrepancy starts to evolve at 50μs approximately.The discrepancy is likely caused by the re f l ections off the components of the setup in the proximity of the emitter and receiver.Since our main concern was to properly reproduce the f i rst arrival and few cycles after,the presence of such deviation was considered within the error of the calibration.

Although both P-and S-wave transducers were used in the laboratory experiments,the f i nal calibration of the numerical model employed S-wave experimental traces only.We treated P-wave calibration results as of limited reliability due to the fact that P-wave transducers were individually produced in the laboratory and the transfer function was not determined.However,we considered this problem less signi f i cant compared to the inaccuracy caused by the presence of one-dimensional(1D)longitudinal stress waves(“bar waves”);see Appendix B for details.Calibration can be improved by utilizing both modes.The resultant best-f i t KPand KScould be then estimated by overlapping GOF plots evaluated for each mode individually.

The low values of attenuation coef f i cients in materials used in calibration(Boyer and Gall,1984;Yang and Turner,2004;Sun et al.,2011)allowed us to neglect the intrinsic attenuation effects throughout the calibration and simulate the fully elastic wave propagation in the calibration setup.Having obtained the values of VPand VS,we performed numerical simulation of elastic wave propagation in a wet bentonite sample.A considerable mismatch was observed in the f i rst 2 cycles(0-55μs)between experimental traces,corresponding to quasi-elastic simulation of S-wave propagation(Fig.5c).The reason for the mismatch is likely the pronounced inelastic behavior of bentonite and resultant attenuation of the signal.It has been noted that clay content in rocks strongly affects attenuation of seismic waves(Barton,2006 and references therein).Literature survey revealed that quality factor of bentonite is three orders of magnitude lower than that of the materials used in the setup(QP,bent≈5,QP,setup≈1000).The simulation of the laboratory experiments in a distinctly viscoelastic material using elastic FD code is likely to yield erroneous results.Thus,to optimize the accuracy of wave propagation in viscoelastic samples, one needs to introduce an appropriate rheological model(e.g. Fig.4b)for the chosen material at certain conditions.However, complex viscoelastic models greatly increase the computational time.Emmerich and Korn(1987)showed that a second-or thirdorder approximation(i.e.a model that includes two or three Maxwell elements)is suf f i cient for most practical applications.The computational time of an elastic subdomain model was approximately 60 s,whereas the viscoelastic simulations(Fig.4)took about 130 s.

The estimation of attenuation was carried out following an iterative process,resulting in a range of GOF minima(red rectangle, see Fig.5).The inconsistent variation of GOF at QP=2 with QSis related to the nature of GOF calculation and may be caused by sudden appearance of spurious peaks in the numerical traces.

The variation of QScauses greater change in the GOF than equivalent QPvariation,as it can be noticed by a more pronounced vertical color gradient compared to the horizontal one(Fig.5a).The actual minimum of the GOF function is smeared over a certain range of(QP,QS)due to limited resolution and the error of the calibration method.However,the more pronounced minimum of the GOF function along a chosen direction should impose the stronger constraint on its exact location in that direction.Consequently,the value of QPobtained through the iterative procedure which employed S-wave traces seems to be of limited reliability with respect to QSand can be estimated only approximately. Therefore,to accurately evaluate QPand QS,both P-and S-wave simulations need to be analyzed individually.The range of possible values of QPand QSthen would be minima of two overlapped GOF plots.Nevertheless,our f i ndings of QP≈10,QS≈6 conform to the data reported by Karakurt(2005)and Astbury and Moore(1970) (QPand QS≈5).Due to resultant high attenuation,there is a large delay(～10μs)between a theoretically predicted S-wave arrival and that predicted by a viscoelastic model when comparing the locations of red triangle and f i rst negative peak(Fig.5b-d).

The work f l ow diagram of the calibration is illustrated in Fig.6. The set of parameters to be varied is the only difference between the calibration and attenuation estimation work f l ows.The main advantage of the iterative routine is that it can be applied to anumerical model regardless of wave propagation solver’s nature (e.g.f i nite difference,f i nite element)and algorithms used in its code.Other advantages of the method are the relative simplicity of the implementation,the computational speed,and the accuracy of the short-term prediction.We focused on the study of the bentonite since the ultimate task of the research is having a tool to simulate wave propagation in engineered barriers.However,the aforementioned work f l ows do not rely on the particular laboratory setup and the material.Therefore,the analogous procedure may be applied to ultrasonic facilities of various con f i gurations,provided that an appropriate mesh is created,following constraints imposed by the chosen numerical method.Moreover,it may be employed to estimate the attenuation of a variety of materials that exhibit viscoelastic behavior(porous rocks,clays,organic materials).The method shows good reliability in resultant values of both VPand VS, although datasets corresponding to only one mode(S-wave)were used.As a result,we have a universal step-by-step guide on the calibration of a numerical solver in order to turn it into a numerical laboratory,capable of reproducing real ultrasonic datasets.

Fig.6.The roadmap of the calibration procedure.

5.Conclusions

The study tries to propose a work f l ow for numerical model calibration,which,if followed,leads to developing a numerical“laboratory”,capable of reproducing real ultrasonic datasets.The aim was also to estimate the attenuation of montmorillonite using the calibrated model.To evaluate the attenuation in a viscoelastic material,the contribution of non-calibrated numerical model to the mismatch should be f i rst eliminated when varying the rheological model and minimizing the mismatch between laboratory and numerical traces.Due to the imperfect contact interfaces in reality,the wave propagation velocities in the material of interest are overestimated.Our investigation revealed that in some cases,the overestimation may be as high as 30%.Therefore,a proper calibration of the numerical tools is of paramount importance for further reliable and realistic simulations.

The resultant calibrated model was successfully applied,allowing the estimation of QPand QSof a 52%water-saturated bentonite. The yielded values con f i rm to those reported in the literature. However the reliability in QPestimation is questionable due to the“available data”limitation.It is shown that strong attenuation properties(QP≈10,QS≈6)of the bentonite must be taken into account in geophysical modeling sensitive to inelastic effects.

Although we were able to consistently reproduce S-wave propagation in standard laboratory tests,signi f i cant discrepancies between numerical and laboratory P-wave traces remained after calibration.This fact indicates that 2D numerical models are likely to oversimplify and disregard some important features of the real wave propagation.Therefore,for higher accuracy and better resemblance,3D models and GSLS-supporting solvers of inevitable high computation cost should be employed.

The main types of host rocks used for the HLRW disposal are granite,clay,salt,and basalt.Therefore,the calibration technique needs to be further tested and validated on the above-mentioned rocks to enable the numerical simulation of monitoring procedure in real geological repository.

Acknowledgments

The authors would like to thank Simon Harvey and Qi Zhao for valuable discussion and advice.This work has been supported by the Natural Sciences and Engineering Research Council of Canada through Discovery Grant 341275(G.Grasselli).

Appendix A.Goodness-of-Fit(GOF)function estimation

The GOF was de f i ned as a measure of the similarity between the numerical and laboratory traces.A Matlab script was written to estimate the mismatch for a set of numerical traces corresponding to a desired laboratory dataset.The user has to specify the time range of interest and the number of consecutive waveform peaks (characteristic points)that lie within the chosen range on a laboratory trace.As we were concerned about the matching in the f i rstarrivals,the set of characteristic points always contained at least one point,corresponding to the f i rst-arrivals registered by the transducer.The user is then offered to manually pick the locations of the points on a laboratory ultrasonic trace(Fig.A1).The chosen peaks ideally should be well-distinguished and correspond to the dominant frequency of the signal.The GOF function for a numerical trace is calculated following the equation:

where Ti，exp（i=1，2，…，n）is the times corresponding to n consecutive peaks of the experimental trace provided by user; Ti，num（i=1，2，…，n）is the times of the consecutive numerical waveform peaks reproducing those picked on the experimental trace.Ti，numis determined by analyzing the numerical waveform either automatically(best for iterative routine over a large number of models)or manually(best for small amount of numerical models).It may be observed that a global minimum of the GOF function(GOF=0)is achieved in the case of perfect match of f i rst n peaks,whereas the maximum of GOF is unconstrained and corresponds to the largest mismatch.This calculation is subsequently applied to all numerical traces related to a chosen dataset and the color plot representing the values of GOF in the space of varied parameters(KP,KSor QP,QS)space is obtained(see Figs.2 and 4).

Fig.A1.Illustration of GOF evaluation:user-de f i ned reference points and corresponding numerical peaks are shown by blue squares and red circles,respectively.The projection of the former and the latter on time axis results in sets of Ti,expand Ti,numrespectively.

Appendix B.Bar waves

The length-to-diameter ratio of the employed laboratory setup was greater than 5.Tisato and Marelli(2013)indicated that bar wave identi f i cation and exclusion were of a primary concern, although SNR was low enough to clearly recognize P-wave arrival. On the contrary,in the case of shear waves,the amplitudes were in most cases signi f i cantly higher than those of bar waves.

The limited ability of 2D numerical model to properly reproduce 3D wave propagation hampers the interpretation of the ultrasonic traces.In the laboratory framework,the signal is re f l ected off the cylindrical boundary.As a result,the conversion of the total input energy into bar waves might be highly underestimated in 2D numerical simulations.The bar waves will affect the traces to an extent as long as the computational domain includes the laterally distributed high impedance contrasts(air-metal bar contrast).

As an attempt to illustrate and justify the hypothesis,we carried out the simulation in the equivalent layered model for both P-and S-wave signals.This model,representing the inner part of the subdomain used above(Fig.1),is unbounded(no air)and layered identically.In this simulation,the bar waves are absent as the medium is transversely isotropic and the air-sample assembly boundary is not presented.The results are compared to those obtained through propagation of S-and P-wave signals in the subdomain model(Fig.A2a and b,respectively).Expected absence of the bar wave effect on S-wave numerical traces is evident:two traces(Fig.A2a)are overlapped within the time range of interest (time before S-wave arrival).Therefore,one should not expect the appearance of strong motion after bar-wave arrival in experimental S-wave ultrasonic traces.On the contrary,there is a large discrepancy between the traces corresponding to bounded and unbounded models in P-wave simulations(Fig.A2b)between P-and S-wave arrivals(20-40μs).Finally,we ran P-wave simulations in a“no-sample”model using KS=0.76,Kp=0.9 and compared the results with the experimental trace(Fig.A2c).The match in f i rstarrivals and one successive cycle(0-35μs)is within the error, whereas the large motion corresponding to the bar wave arrival seems to be absolutely missing in the numerical trace.Thesephenomena lend validity to the aforementioned hypothesis that 2D code is likely incapable of proper numerical depiction of the refl ections off the cylindrical boundary that takes place in reality and underestimates the energy carried by the bar wave.

Fig.A2.Numerical traces corresponding to S-and P-wave experiments in bounded(a)/unbounded(b)models and“no-sample”model(c)for P-wave propagation using KP,KSobtained from calibration(top to bottom,respectively).Red circles,green triangles and blue squares show analytically predicted P-,bar,and S-wave arrivals,respectively.The mismatch after two cycles(bottom)is likely to be related to the“bar wave”effects.

Appendix C.Errors of the calibration

It is not straightforward to list all factors that may contribute to the error of both the calibration and attenuation estimation.To the aforementioned errors related to resolution in KP,KS(or QP,QS)and smear of GOF minimum location,one should add the uncertainty caused by the thickness and coupling of transducers.Piezoceramic discs utilized in P-and S-transducers were 3 mm thick(22 grid points,7%-10%of the total sample assembly length),which imposes uncertainty to excitation and acquisition location.In the numerical model,the signal was excited at a horizontal set of nodes corresponding to the transducer surface that is in contact with the assembly.In addition to that,the numerical input and output were assumed to be non-coupled and perfectly 1D(excitation in vertical or horizontal direction for P-and S-waves,respectively).In reality, the disc is subjected to a complex,often coupled deformation. However,these issues were partially taken care of during the calibration and should be of second-order effect on the f i nal results of attenuation estimation.

The oscillations between 0 and 15μs(Fig.5b-d)are caused by electromagnetic disturbances generated by the emitter and captured instantaneously by the receiver and are not related to wave propagation in montmorillonite.

Con f l ict of interest

We wish to con fi rm that there are no known con fl icts of interest associated with this publication and there has been no signi fi cant fi nancial support for this work that could have in fl uenced its outcome.

Alexander W,McKinley L.Deep geological disposal of radioactive waste(radioactivity in the environment).Rotterdam:Elsevier;2007.

Alford RM,Kelly KR,Boore D.Accuracy of f i nite-difference modeling of the acoustic wave equation.Geophysics 1974;39(6):834-42.

Alonso EE,Springman SM,Ng CWW.Monitoring large-scale tests for nuclear waste disposal.Geotechnical and Geological Engineering 2008;26(6):817-26.

Astbury NF,Moore F.Torsional hysteresis in plastic clay.Rheologica Acta 1970;9(1): 124-9.

Auerkari P.Mechanical and physical properties of engineering alumina ceramics. Espoo,Finland:VTT Manufacturing Technology;1996.

Barton N.Rock quality,seismic velocity,attenuation and anisotropy.1st ed.New York:CRC Press;2006.

Birch F.The velocity of compressional waves in rocks to 10 kilobars:1.Journal of Geophysical Research 1960;65(4):1083-102.

Bohlen T.Parallel 3-D viscoelastic f i nite-difference seismic modeling.Computers and Geosciences 2002;28:887-99.

Boyer HE,Gall TL.Metals handbook:desk edition.Ohio:Asm.Intl.Metals Park; 1984.

Chapman N,McCombie C.Principles and standards for the disposal of long-lived radioactive wastes.Rotterdam:Elsevier;2003.

Emmerich H,Korn M.Incorporation of attenuation into time-domain computations of seismic wave fi elds.Geophysics 1987;52(9):1252-64.

Karakurt N.Estimating attenuation properties of bentonite layer in cut bank oil fi eld,Glacier county,Montana.Texas:Texas A&M University;2005.PhD Thesis.

Lajudie A,Raynal J,Petit J.Clay-based materials for engineered barriers:a review. In:Scienti fi c basis for nuclear waste management(Part 1),materials research society symposium proceedings.Pittsburgh,USA:Materials Research Society; 1994.p.221-9.

Lakes RS.Viscoelastic materials.London:Cambridge University Press;2009.

Liu HP,Anderson DL,Kanamori H.Velocity dispersion due to anelasticity:implications for seismology and mantle composition.Geophysical Journal International 1976;47(1):41-58.

Madonna C,Almqvist BSG,Saenger EH.Digital rock physics:numerical prediction of pressure-dependent ultrasonic velocities using micro-CT imaging.Geophysical Journal International 2012;189(3):1475-82.

Manukyan E,Maurer H,Marelli S,Greenhalgh SA,Green AG.Receiver-coupling effects in seismic waveform inversions.Geophysics 2012;77(1):R57-63.

Marelli S,Manukyan E,Maurer H,Greenhalgh SA,Green AG.Appraisal of waveform repeatability for crosshole and hole-to-tunnel seismic monitoring of radioactive waste repositories.Geophysics 2010;75(5):Q21.

O’Connell RJ,Budiansky B.Measures of dissipation in viscoelastic media. Geophysical Research Letters 1978;5(1):5-8.

Saenger EH,Madonna C,Frehner M,Almqvist BSG.Numerical support of laboratory experiments:attenuation and velocity estimations.Acta Geophysica 2014;62(1):1-11.

Sun EW,Cao WW,Han P.Frequency dispersion of longitudinal ultrasonic velocity and attenuation in[001]c-Poled 0.24Pb(In1/2Nb1/2)O3-0.45Pb(Mg1/3Nb2/3)O3-0.31PbTiO3single crystal.IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control 2011;58(8):1669-73.

Tisato N,Marelli S.Laboratory measurements of the longitudinal and transverse wave velocities of compacted bentonite as a function of water content,temperature,and con fi ning pressure.Journal of Geophysical Research:Solid Earth 2013;118(7):3380-93.

Villar MV,Martin PL,Barcala JM.Modi fi cation of physical,mechanical and hydraulic properties of bentonite by thermo-hydraulic gradients.Engineering Geology 2005;81(3):284-97.

Yang L,Turner JA.Attenuation of ultrasonic waves in rolled metals.The Journal of the Acoustical Society of America 2004;116(6):3319-27.

Anton Biryukov is enrolled in Master of Applied Science in Civil Engineering program at the University of Toronto. He obtained Bachelor of Applied Physics and Mathematics degree from Moscow Institute of Physics and Technology. He is currently working as a research and teaching assistant at the Department of Civil Engineering.His research interests cover experimental subsurface geophysics, microseismicity induced by oil and gas exploration and reservoir stimulation,numerical geophysical and geomechanical modeling using f i nite-difference/combined FEM-DEM tools.He has participated in several geophysical and rock engineering conferences,presenting studies on microseismic source location and numerical modeling of laboratory experiments.Since 2013,he is an active member of several academic organizations,including Canadian Society of Exploration Geophysicists(CSEG),Canadian Society of Petroleum Geologists(CSPG),Civil Engineering Graduate Students Association and University of Toronto Education Workers.

*Corresponding author.Tel.:+1 647 771 2953.

E-mail address:anton.biryukov@mail.utoronto.ca(A.Biryukov).

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

1674-7755?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

http://dx.doi.org/10.1016/j.jrmge.2014.10.002