Excavation Damage Zone fracture modelling for seismic tomography: A comparison of explicit fractures and effective medium approaches

We model the full wavefield produced by a seismic velocity survey and optimise the representation of the fracture zone to best match field waveforms. The velocity survey was part of a mapping study on fractures in the Excavation Damage Zone (EDZ) of ONKALO underground research facility at Olkiluoto. The EDZ results from excavation of the rock mass, which modifies stress conditions changing the nature and behaviour of pre-existing fractures and generating new fracturing. These fractures act as the main transport pathways for contaminants both in and out of a geological disposal facility (GDF). Our goal is to test different representations of the fracture zone and to determine which models most successfully improve the interpretation of the fracture zone, producing estimates of a key unknown parameter, fracture stiffness, in addition to fracture sizes, fracture geometry, fracture density and crack density. We use modelling techniques previously tested in theoretical and laboratory studies and assess their performance on a real engineering problem. The paper introduces the field experiment and relevant information from the GDF in Finland. It describes the methodologies used for representing the fracture networks in the models — Explicit Fracture models with two approximations called Pixelised Fracture Model (PFM) and Equivalent Discrete Fracture Medium (EDFM), the Effective Medium (EM) model, and two versions of the Localised Effective Medium (LEM) model (LEM fine, LEM thick). These alternative representations were used within models of the field experiment and the calculated waveforms were used in an iterative inversion for fracture stiffness. Results show that the EM model and the EDFM model were unsuccessful in matching recorded waveforms. The fine LEM model and the explicit PFM model produced the best results especially after iterative optimisation of the fracture stiffness, giving confidence that further optimisation will lead to improved characterisation of the fracturing from the full waveform data.


Introduction
Identifying and characterising the properties of fractures and fracture zones is of widespread importance in a number of rock engineering applications both for assessing stability and for identifying fluid pathways (e.g., Refs.1-4).Such fracture zones directly modify seismic waves and seismic velocity scans are therefore used to provide data on depths and characteristics of fracture zones (e.g., Refs.5-8).There is ambiguity in such interpretations on the density, sizes, orientations and stiffnesses of fractures that have caused the bulk changes in the seismic waveforms (e.g., Refs.9-13).We seek to establish how best to use complexity of such fracture zones which challenges the desired combination of efficiency and accuracy in the fracture representations.To assess the performance of different fracture representations on fieldscale data, we selected a velocity scan experiment performed as part of a mapping study on fractures in the Excavation Damage Zone (EDZ) of ONKALO underground research facility at Olkiluoto, the future GDF for nuclear waste in Finland.Nuclear waste must be safely disposed of in a permanent geological disposal facility (GDF) when the temperature of the spent fuel drops to safe levels.At present, countries producing highlevel waste (HLW) store the spent fuel in cooling water pools, while researching and designing a permanent GDF solution.Posiva Oy is a company responsible for the design and construction of the ONKALO underground research facility (ONKALO means ''cavity'' in Finnish).In 2007, Posiva Oy carried out an EDZ research programme aimed at best characterising the EDZ.As part of this, a cross-hole velocity tomography study was carried out by ASC Ltd. 8,16 to measure seismic wave velocities, in the testing tunnel of ONKALO ONK-TKU-3620 niche at 345 m depth (see Fig. 1).The survey concentrated only on seismic velocities not the full waveforms.Variations in the seismic wave velocities were used to interpret damage in the rock.In this work, we create models of the fracture zone and model the velocity surveys from the above study using the recorded waveforms for comparison with those produced by the models.Schoenberg 5 introduced the linear slip model based on displacement discontinuity conditions at an interface as a model of a fracture and demonstrated its frequency-dependent effects on seismic waves.Pyrak-Nolte et al 17 demonstrated through laboratory experiments that the displacement discontinuity model is appropriate for seismic wave transmission and reflection from natural fractures.Hildyard 18 demonstrated that an explicit model of a fracture with two fracture surfaces implementing displacement discontinuity conditions was successful in matching waveforms from both laboratory and in situ recorded waveforms.Parastatidis et al 14 investigated three different representations of fracturing and assessed their performance in matching waveforms from a laboratory experiment.He used an Effective Medium (EM) model, an explicit fracture model, and a Localised Effective Medium (LEM) model.He confirmed conclusions from Hildyard 18 that waveforms from the explicit fracture models matched those from the experiment and found that the Localised Effective Medium (LEM) model produced equally good results to the Explicit Fracture Model (EFM) provided the LEM layer was sufficiently thin.The Effective Medium (EM) model was unsuccessful at matching the experiment unless all results were filtered to very low frequencies.As the LEM layers were made thicker their accurate frequency range reduced until results matched the EM model.This work investigates the performance of these three types of fracture zone representations.However, the increased geometrical complexity of the in situ fracturing means that different approximated versions of explicit fracture models were developed along with different variations of the LEM model leading to a wider range of models.

The experiment
Since 2004, Posiva Oy has carried out tests and experiments in the future host rock in Olkiluoto, to have a better understanding of the geology and mechanical behaviour of the rock in this area.In 2015, the company obtained a licence to move forward from the testing stage to the construction of the GDF.Further important research is to understand the development of fractures during the excavation process when constructing the GDF and the fracture behaviour during the closing of the tunnels in response to the back-fill material and the thermal loading from the spent fuel.The excavation process or any future changes in the bedrock stresses due to external processes can result in irreversible changes developing in the rock mass such as rock damage or in initiation of new fractures and deformation of the existing fracture network.The regions directly surrounding man-made openings are referred to as the Excavation Damage Zone (EDZ) and the Excavation disturbed Zone (EdZ) (Fig. 2).Tsang et al 1 reviewed the hydromechanical processes linked with man-made openings and provided a definition and details of the evolution of the EDZ and the EdZ for all stages from the construction of the GDF to its longterm safety after closure of the GDF.The conclusions from Ref. 1 on crystalline rocks were based on experimental and modelling studies from several Underground Research Laboratories (URL).The EdZ is a zone some distance from the tunnel surface where the changes are stress related and are reversible having no negative effect on the GDF safety for waste isolation.For crystalline rock it is difficult to define the outer boundary of this zone. 1 The EDZ for crystalline rock is a region with irreversible changes and 'freshly opened' fractures.The thickness of the EDZ depends on the tunnel excavation method.Siren et al 2 has grouped the EDZ based on the origin of the damage into two subcategories.The first links the origin of the fracture opening to the excavation phase itself and is called construction-induced EDZ (EDZ  ) and these fractures are closer to the surface of the tunnel (Fig. 2).The second group is related to stress changes.It is called stress-induced EDZ (EDZ  ) and it is further from the tunnel surface (Fig. 2).The EDZ is stable for crystalline rock after the construction stage (of the GDF). 1 During the early closing stage, water from the rock will resaturate the EDZ, while the temperature from the spent fuel will have a dynamic impact on the behaviour of the EDZ.The EDZ fractures will begin to close as the backfilling material swells and applies pressure.Tsang et al 1 concludes that the relatively high permeability of the EDZ must not be studied separately but as a total flow system, since when the high permeability zone is surrounded by low permeability, the water supply in the EDZ will be poor.
In 2007, Posiva Oy carried out an EDZ research programme to study how the excavation process itself affects the effectiveness of the rock by creating new possible paths for fluid flow.The EDZ experiment was conducted in the ONKALO underground research facility at Olkiluoto, before it was expanded into a GDF (Fig. 1).The key objectives of this programme were to create a method for characterising the EDZ, to improve the excavation methods to minimise damage, to evaluate flow paths for nuclear waste leakage and to identify solutions to eliminate migration (e.g., Refs.2-4, 7, 8).Several monitoring methods have been evaluated to monitor the evolution of an EDZ such as microseismic and acoustic emission (AE) monitoring (e.g., Refs.9-11), and cross-hole tomography (e.g., Refs.12, 13).As part of this EDZ programme, ASC Ltd. 8,16 carried out a cross-hole velocity tomography study to measure seismic wave velocities in the testing tunnel of ONKALO ONK-TKU-3620 niche.The variations in the seismic wave velocities were used to interpret the damage in the rock due to fractures which might have been created during excavation near the tunnel surface.The ONK-TKU3620 was excavated in exceptionally good rock conditions, in the absence of intense natural fracturing or intersections of brittle deformation zones. 7,16,20The survey focused on seismic velocities obtained from first arrivals not on examining the full waveforms.However, the modelling study in this paper makes use of the full waveforms from this velocity survey.Several reports produced by Posiva Oy include a detailed description of the geology of the surrounding area of the testing tunnel (Fig. 3(a)).These studies used a combination of methods to map the lithology, in-situ stress and fracture orientation including well logging, geophysical measurements, and photogrammetry.

Background lithology of the studied area
The main rock types in the tunnel region are veined gneiss (VGN) and pegmatitic granite (PGR) (Fig. 3(a)).Fig. 3(a) shows the mapped lithology of the testing tunnel ONK-TKU-3620 niche and Fig. 3(b) shows a plan view of the studied area for the P-wave velocity survey.The rock type in the survey area is VGN (Fig. 3(a)).The average VGN dip and dip direction is 45∕163 • . 7,19The image shows the extent of ONKALO that was available during the execution of the EDZ experiment, before being expanded into GDF.

3D mapped fractures and orientation
There are three sets of fractures which pre-date the excavation.The orientation (Dip and Dip Direction) of these fracture sets is 38∕159 • , 85∕267 • and 84∕344 • .The area where the seismic tomography has taken place has a detailed fracture model which has been created by combining different geophysical methods, well logging and sawing of blocks E4 to E7 as shown in Fig. 3(b).Using the results from the sawing blocks, it was possible to identify the fractures in the other geophysical data.
In the current fracture model, the fractures are classified according to their origin with respect to the natural fractures, Natural fractures are opened by excavation, EDZ  fractures, and one fracture F1 classified as EDZ  . 20The fracture model in the sawed area is presented in Fig. 3(c). 7

Survey setup
The P-wave velocities measured during the active velocity survey were acquired using ASC's ultrasonic and acoustic emission (AE) monitoring equipment system.This consisted of an array of twenty-four  this surface array were coupled to the ground surface using rock bolts.
For each single active velocity survey, one transducer acts as an active source where it is pulsed with a high voltage excitation signal (provided by the ASC pulser system).The response of the excited transducer propagates a signal through the surrounding rock mass which is detected by the remaining AE transducers in both the borehole and surface sensor arrays.The amplified AE transducer signals are measured by the data acquisition system which records the full waveform data.Each of the AE transducers in the sensor array in turn acts as an active source.Therefore for each borehole pair there are 24 active shots (24 shots × 23 recordings) with 4 survey instalments taking place on each 2D frame (Fig. 3(d)).This creates a total of 552 acquired waveforms per survey and a total of 2208 waveforms per 2D frame (since there are four surveys per frame as shown in Fig. 3(d)).For the first survey, both frames were placed at a depth between 0 m-0.4 m.In the second survey, the first frame was placed at a depth between 0.4 m-0.8 m with the second frame remaining in the initial position.For the third and subsequent surveys, the depths of both frames were varied and increased as required.The data from the AE transducer surface array had a low signal to noise ratio 8 and was not used in this study for comparisons with the models.

Modelling the experiment
The aim of this work is to study wave propagation in various directions and compare the model with the acquired survey data.Therefore it is necessary to create a detailed and accurate representation of the modelled area.The most essential information is the lithology, the rock properties based on dynamic measurements, the stress state, an accurate fracture model, and a clear seismic source representation.From the reports provided by Posiva Oy and ASC Ltd., 7,8,19,21,22 we have detailed information on the lithology, stress state and fracture model of the studied area.The active seismic source used for each of the surveys carried out must be retrieved by a source inversion technique similar to the one described in Ref. 15.The studied volume is 6 m 3 large, covering blocks E3 to E8 as shown in Fig. 3(b).The complete seismic velocity survey data acquired during the active survey study, consist of tens of thousands of waveforms.The models are computationally expensive and it is impractical to model the full survey volume.Instead we model a region of the survey (1 m 3 ) and compare the model results with a selected subset of the survey data.The region modelled was selected to avoid geometrical complexity to most easily fit orthogonality restrictions for explicit fractures in WAVE3D. 18,23,24AVE3D uses a staggered finite difference grid, which is widely used to solve wave propagation numerically both in elastodynamics and electrodynamics.It solves to fourth-order accuracy which reduces numerical dispersion and increases the usable frequency range for a given element size ultimately increasing efficiency.

Design of the models
This problem needs to be modelled in three dimensions firstly to obtain the correct wave physics and also to study the wave propagation in all the possible azimuths of the survey.In designing the models we need to consider numerical restrictions.The first restriction relates to the wave frequency and the element (or zone) size to avoid numerical dispersion.We constrain this using a conservative relationship based on second-order accuracy ( ≥ 10 x, where  is the wavelength and x is the element size), since not all regions in the models are solved to fourth order accuracy.Based on the frequencies observed in the recorded waveforms, the predominant period of the wave is 11.8 μs and based on the report 8 the minimum velocity is 4800 m/s.This leads to an element size no larger than 5.6 mm and an element size of 5 mm was selected.A volume of 1 m 3 consists of 8,000,000 elements.A further padding of 100 elements was used on each side of the volume making the models significantly larger.The modelled region requires a detailed and verified fracture model which are the blocks labelled E4 to E7 which were sawed so that fractures were mapped.The lithology of the survey region is close to a contact between pegmatite Granite and veined Gneiss.Since our focus is on the effects on wave propagation due to fractures, it is important to keep the lithology as simple as possible.Based on the detailed mapping of the lithology, 7 a region with minimal mixture of different rock types was selected which is block E4 consisting of veined Gneiss (Fig. 3(b)).This block is 1 m 3 in size with four boreholes on the edges of the block (Borehole numbers 34, 35, 43 and 44 in Fig. 3(b)).For each block, there are six complete surveys based on the crosshole borehole pairs (34-43, 34-35, 34-44, 35-43, 35-44 and 43-44).For each of these surveys there are 552 waveforms.The model needs to focus on specific waveforms from particular source-receiver pairs.Note that for convenience in the paper we refer throughout to these source-receiver pairs as ''raypaths'' even though we are modelling the full wave behaviour and not making approximations with ray theory.

Source inversion
The model requires an appropriate repeatable source pulse that best matches the excitation induced in the rock by the transducer.The reports on the survey, 8,21 have detailed descriptions of the data acquisition system and the survey design, but, there is no detailed information describing the active source as used for each survey The only information is provided by the waveform recorded at the transmitting sensor.However, the received signal at the sensor position is more like a square wave due to 'clipping' due to amplitude saturation, making any frequency analysis difficult.The clipping is observed for all surveys for the sensor next to the source transducer.Clipping also occurs for the two sensors above (10 cm upwards) and the two below (10 cm below).Such waveforms are not required in our analysis which focuses on waveforms received at the opposite borehole where the raypaths cross the fractures.It does however limit the useful source information that we can get from these particular waveforms.From these recordings though, we concluded that the function of the effective source is not a single spike but more of a double ''Ricker-like'' shape with at least three full cycles.We therefore need to invert for the source based on waveforms recorded at some distance from the source.An inversion process based on deconvolution of the transfer function has been successfully applied to data from a laboratory experiment. 14owever, in that experiment, a source could readily be inverted from a base case where waveforms were recorded in a solid homogeneous steel block.In this survey though, there is nothing equivalent to this.However, based on the fracture model from the sawing and the mapped fractures from the boreholes, 7 the surface for the block E4 has no fractures and can be considered as homogeneous.We used waveforms from the survey within this relatively homogeneous region and inverted them using the same process as Ref. 14.However, since the rock is not completely isotropic and homogeneous, the source inversion process is more complex than in Ref. 14 and cannot be inverted uniquely from a single shot position.Instead, the source inversion process was repeated several times for different borehole pairs in the block.Seven raypaths were selected within the non-fractured region of the block, all of them at a depth of 0.7 m to 0.8 m.The raypaths are for specific borehole pairs and the inversion was applied in both directions giving fourteen waveforms.These waveforms were found to be similar in arrival time, frequency, and amplitude.They were selected based on having limited high-frequency noise, and to have a horizontal raypath where the receiver is at the same depth as the source transducer in the opposite borehole.For the inversion process, we modelled these using an initial arbitrary source as a 0.8 MHz frequency periodic wave within a homogeneous isotropic medium with dynamic material properties as defined in the report. 21Source waveforms were inverted for each of all the raypaths and these were phase shifted to best match one another.The fourteen sources were then combined to create a mean source.In

Table 2
Maximum cross-correlation coefficient between survey and modelled data with the inverted source and the time lag for the selected raypaths.These raypaths have no presence of fracture and were considered as homogeneous in order to apply the source inversion method.The high value of the cross-correlation coefficient along with the small time lag (in μs) shows that the source inversion method has produced a reliable source.the mean source we applied a tapering function to include only the first part of the source, as the rest might be due to reflections and ringdown count.The mean source was then tested in homogeneous models for the same positions as the waveforms used for the inversion process.

Max. coefficient
The correlation coefficients along with the time shifts needed to reach maximum correlation were calculated (Table 2).Based on these values we conclude that the derived source is sufficiently reliable for use in the full velocity scan models.

Implementation of fracture zone
Posiva Oy provided us with the design files of the detailed fracture models (Fig. 3(c)).Fig. 4(a) is a sketch of a 2D cross-section of block E4 for the borehole pair B34-B43 from the velocity surveys showing the fractures.The surface area of each fracture within the area of interest has been measured separately.The fracture area was then used for the design of the models.However, most of the fracture surfaces are not planar and further information was utilised from the Report 7 based on the stereo-nets of the measured fractures, in order to describe each fracture with only two angles (dip and dip direction).This information is summarised in the Appendix C Table 4.The majority of the fractures have a smooth dipping angle between 6 to 20 • .

Fracture representations: explicit models
There are two main methods for representing fractures in numerical models.The first is to use displacement discontinuities to represent explicit fractures within a background medium.The displacement discontinuity model proposed by Schoenberg 5 considers two elastic homogeneous and isotropic spaces in a non-fixed contact of zero thickness connected by linear springs.Explicit fractures are implemented in WAVE3D as two explicit surfaces with displacement discontinuity conditions coupling the surfaces and with contact conditions at fracture edges. 18The fractures can have linear normal and shear fracture stiffness or non-linear stress-dependent fracture stiffness.The fractures also have the capacity for shear and tensile failure but these properties are not relevant to this study.Full details on the implementation are available in Refs.18, 24.WAVE3D uses an orthogonal staggered grid and this limits its ability to represent complex fracture geometries.However, Parastatidis 15 describes two approximate implementations for a dipping fracture using the explicit fracture model in WAVE3D.
The first implementation is to ''pixelise'' the fracture surface, or the Pixelised Fracture Method (PFM).This approach is simply to ''pixelise'' the dipping fracture with smaller horizontal and vertical fractures.The ratio of the size between horizontal and vertical fractures depends on the angle of the fracture (e.g. if the dip angle is smooth the size of the vertical fractures is smaller than the size of the horizontal, see F1 fracture in Fig. 4(b)).
The second approach is based on the crack density  as described in Appendix A Eq. ( 3), with the creation of random horizontal and vertical fractures with the same total crack density as the single fracture.This approach could be considered as an equivalent discrete fractures medium (EDFM).The equivalent discrete fracture medium (EDFM) is an explicit fracture model which, instead of mapping the exact fracture geometry and orientation, uses orthogonal vertical and horizontal fractures at random positions in the same volume as the original fracture.These random fractures have the same total crack density  as the single fracture with complex geometry.Both methods were validated and guidelines were reached on how to represent the dipping surfaces for a single fracture model. 15Those conclusions are used here to construct the explicit fracture models.There will be two cases of explicit fracture models for each survey; the PFM and the EDFM (Figs. 4(b) and 4(c)).

Fracture representations: Effective medium models
In modelling seismic waves, a common approach for fracture representation is an effective medium model (EM) (e.g., Refs.15, 25-27).The effective medium (EM) combines the effect of the fractures and the rock into a single medium.For example, if a medium contains a single set of parallel fractures, this leads to the simplest effective medium which can be described by five elastic constants (Appendix A, Eq. ( 6)), provided the wavelength is sufficiently large (at least an order of magnitude larger than the size of the fractures).The case of a single set of parallel fractures creates a transversely isotropic medium.The value to describe the fracture population in the EM model is the cracks per unit length parameter, 1/L (Appendix A, Eq. ( 5)).''L'' is the total area of cracks in the volume of interest.The EM stiffness matrix also depends on the orientations of fractures, the normal and shear fracture stiffness, and the background material properties.For the EM model, the total area of the fractures within the total volume was used to calculate 1/L, and the average dip and dip direction angles was used. 7s a result, the EM model was quite simple to design (Fig. 4(d)).
A further approach based on Effective Medium theory is the Localised Effective Medium (LEM).The LEM model uses the same approach to constructing the stiffness matrix as the EM model (Appendix A, Eq. ( 5)).The major difference between EM and LEM is that in the case of EM the fractures in the stiffness matrix are represented as a sum of all the fracture areas in the modelled volume (see Appendix A Eq. ( 5)).In contrast, the LEM solves the five constant effective medium moduli locally for material surrounding each fracture.The material properties (EM stiffness matrix) are calculated separately for each zone containing part of a fracture and hence these properties can vary zone by zone (e.g., Refs.15, 28-31).This approach can be seen as a hybrid between the explicit fractures approach and the Effective Medium approach.Parastatidis et al 31 studied the performance of the LEM approach and concludes that when the LEM layer is very fine the resulting waveforms are close to the explicit model while the thicker the LEM layer is the closer results are to the EM model.
As for the explicit models, two cases of the LEM model were explored.The first case (LEM fine layer) uses a cuboid-shaped volume of LEM material around the fracture with coordinates such that the fracture fits within this volume.The volume and the surface area for each fracture gives a different 1/L value for each fracture region.For each fracture, the fracture angles are also used to rotate the stiffness matrix for the correct orientation of the fracture (Fig. 4(e)) again leading to different stiffness matrices for each fracture region.In the second LEM case (LEM thick layer), the input values are the four corner coordinates of a rectangular fracture and the two angles, then WAVE3D creates a thin LEM layer with a thickness equal to one element while 1/L is equal to 1/x (Fig. 4(f)).The first case produces lower 1/L values compared to the second case.As shown in Ref. 31 the lower the value of 1/L the closer the model is to the EM model.The first case of LEM model still differs from the EM model because, even though the 1/L will be low, the dip and the dip direction angles will be different for each fracture and thus each fracture region has a different stiffness matrix compared to the EM model which is homogeneous.The fracture normal stiffness K  used in the initial model is 3.3 × 10 10 Pa/m (Appendix C Table 4).The values of fracture stiffness were estimated and scaled  according to the size of the fracture as described by Worthington and Lubbe. 32The shear stiffness K  was set at one third of the normal stiffness.These values were selected for the initial models.Based on the results, the stiffness of some fractures were adjusted and rerun in an iterative optimisation in order to better match the data.These five fracture models, the two explicit fracture models (PFM and EDFM), the EM model, and the two LEM models (LEM fine, LEM thick), were inserted in the model of the seismic tomography.,The performance of the initial models will be assessed against the survey data in the next section.

Preliminary model evaluation
Overall, modelling the survey poses many challenges.First of all, there could be a background anisotropy in wave velocities since the VGN formation has a dip and dip direction creating a foliation which might affect the velocities.The source function creates uncertainty as to whether the source from the inversion process is representative in all cases.The source issue is linked to a further problem; in that some sensors of the array may have a better coupling with the surface of the borehole creating differences in amplitude.The coupling between the borehole surface and the receiver and/or transducer have an effect on the transfer function in the acquisition system 13 and can result in differences in amplitude from 5 to 25%.For example in the 'B34-B43 0.8 m' case the recording on borehole 1 (red) has 25% higher amplitude compared to the recording on borehole 2 (blue), even though both waveforms are for propagation along the same path in different directions and hence expected to be identical.The presence of water in the boreholes 8,16 creates further complexity with expected attenuation and slowing of wave velocities.In designing the models, assumptions needed to be made about the geometry of the fractures for the different fracture representations as explained in previous sections.Finally, although fractures were mapped, the fracture stiffness is an important but unknown parameter and needs adjustment for waveforms to match the survey data.Due to the above complexity and to minimise the amount of data, we concentrate on specific raypaths which we will try to optimise in the second stage to improve the model outputs.We therefore focus on waveforms between boreholes 34 and 43.
Between boreholes 34 and 43 there are two main fracture sets.The first one is the F1 fracture which is cutting the whole block with 6 • average dipping and 301 • dip direction.The minimum and maximum depth of the fracture is at 46 cm and 70 cm respectively.The initial fracture stiffness used is 3.33 × 10 10 Pa/m (Appendix C Table 4) and it is based on the size of the fracture which is greater than 1 m 2 .The F1 fractures generated due to the fact that the prevailing stress field was re-distributed around the tunnel opening after the excavation of the void space.This stress re-distribution sometimes exceeds mainly the tensile strength of the rock at the tunnel floor, thereby creating large EDZ  fractures which were studied in detail during the testing in Olkiluoto. 20The second fracture set has been characterised as EDZ fractures.It has been split into four fractures, as shown in Appendix C, Table 4, fractures 9,11,12 and 15.The dip of the fractures varies from 6 • to 12 • and dip direction is approximately 180 • .The size of the fractures is between 0.5 m 2 to 0.017 m 2 with fracture stiffness 5.71 × 10 10 Pa/m to 5 × 10 12 Pa/m.The dip angle of all the fractures presented between boreholes 34 and 43 are almost horizontal.The two fracture sets are highlighted in Fig. 4(a).
From the initial results we selected a specific raypath to demonstrate how comparable models were with the survey data.The selected raypath is the shot on borehole 43 at 0.60 m depth and recording at depth 0.55 m (further examples with different raypaths are presented in Appendix B).The waveform in this case propagates with 3 • angle from horizontal.The raypath crosses the F1 fracture.
Fig. 5 shows waveforms for the five model types for the raypath (Fig. 5(a)) which crosses the F1 fracture.For this raypath, based on travel time calculations, the part of the waveform after the first 0.25 ms is the result of a strong reflection from the EDZ fractures above the F1 fracture.
The first visual impression from the results of the initial models is that the explicit PFM model (Fig. 5(b)) and the thick LEM model (Fig. 5(e)) perform in a similar way and are very close to the survey data.The thick LEM model in this cases works better than the explicit model (Fig. 5(a)).The fine LEM model (Fig. 5(f)) with the high 1/L value in this case causes more attenuation compared to the thick LEM model which is closer to the recording.The EM model (Fig. 5(d)) and EMDF model (Fig. 5(c)) are not at all comparable with the survey data and can be eliminated from further analysis.

Optimisation
Apart from the cases shown above where one or more waveforms from the different fracture models match the survey, the majority suffer from either earlier or later arrival and higher or lower amplitude.Increasing or decreasing the fracture stiffness could have an effect on the final result.
Based on the comparisons of the models against the survey data, we conclude on which fractures need adjustment in stiffness.For the survey in boreholes 34 and 43 (Fig. 5(a)), the waves are more attenuated in terms of amplitude for the paths going through the F1 fracture.The F1 fracture has some unique characteristics -it is large enough to cover the whole block and is smoothly dipping.The depth of the fracture is between 0.5-0.7 m and it is isolated from the other fractures and the free surface.The fact that it is isolated is helpful at this stage since the stiffness of a single fracture will need to be adjusted avoiding any complex system with more fractures which might create non-unique results.

Stress-dependent stiffness versus manual iterative optimisation
Fracture stiffness depends on stress state and the normal fracture stiffness scales with the stress normal to the fracture surface. 33One approach to optimise fracture stiffness is to solve the stress state and to directly link fracture stiffness to the localised stress state.This can lead to a better match and a better physical explanation of the measured waveforms. 15,18In the above studies there was a significant variation in the stress field across the length of the fractures but this is not expected in our case.Since the fractures on the borehole frame 34-43 are almost horizontal (Appendix C Table 4 for fracture dip), the normal stress is close to the vertical stress, and near the tunnel surface this stress will be close to zero.As a result, the stress variation will only cause a small variation in stiffness along the fracture.Hence for this first analysis we instead used manual iterative optimisation to find a single value of stiffness for the entire fracture rather than accounting for variation along the fractures.Such an approach is also an order of magnitude less computationally demanding as the models do not need to also solve the stress state, and multiple cases can be run in parallel.Nevertheless, a stress dependent approach may be an interesting extension in future studies.

Manual iterative optimisation
Based on the observed scaling between fracture size and fracture stiffness 32 the preliminary models were set to use the lowest stiffness values in the observed range.For the size of F1 fractures, the suggested stiffness can range up to a hundred times higher.As a result, we created thirteen cases with varying stiffness (as summarised Fig. 6).For the first ten cases, we increase the normal stiffness from ten times (case 1) to one hundred times (case 10) the initial stiffness value at intervals of 3.33 × 10 11 Pa/m while the shear stiffness is set as a factor of 0.33 of the normal stiffness.For fracture F1,   = 3.33 × 10 10 Pa/m in the initial model, and varies here from   = 3.33 × 10 11 Pa/m to 3.33 × 10 12 Pa/m.In the last three cases, the normal fracture stiffness is constant and the factor for the shear stiffness varies as 0.5, 0.6 and 0.2.We selected a specific raypath to optimise the fracture stiffness.The raypath goes from borehole 43 to 34 with the shot position at 0.6 m depth and recording at 0.55 m on borehole 34, see Fig. 5(a).The raypath only crosses a single fracture (F1) and consequently the waveform is not affected significantly by any other fractures which would complicate the optimisation process.Furthermore, both the fracture and the raypath are relatively deep and away from other fractures near the surface that could add complexity to the waveform, through reflections and diffractions.Fracture F1 is a large fracture that cuts through the whole block at a smooth angle and, as a result, the chosen raypath propagates with an angle of approximately 10 • to the fracture.This path showed significant changes in the waveform for different fracture stiffnesses.The optimisation process was applied to the explicit PFM model and both the fine and course LEM models.A large number of waveforms were produced during this optimisation.We present figures from just four of the thirteen cases.The results of the three models for the chosen raypath are presented in Fig. 8.For this frequency and for these propagation angles relative to the fracture, the small values of stiffness in the first cases, were not enough to change the waveforms for the explicit and for the fine LEM models (Figs.8(a), 8(b) and 8(c)).As the stiffness reaches forty times the initial stiffness, the amplitude of the model waveform is similar to the survey data for a time window between 0.15-0.25 ms (Figs.8(d), 8(e) and 8(f)).When the stiffness goes above fifty times the initial stiffness, the amplitude of the models is much higher than the data for all three models (Figs.8(g), 8(h) and 8(i)).The velocity (from arrivals in the waveforms) does not change significantly for the two LEM cases with maximum differences of up to 2.2 μs.For the explicit case, the velocity remains the same for all different situations.Cases eleven to thirteen, where only the shear stiffness changes, have lower impact on the amplitude of the waveform compared to the changes due to normal stiffness.In this case, the amplitude of the waveform is reduced by about 5% compared to the case five with the same normal stiffness.

Results -Initial vs. ''Best iterative'' cases.
Fig. 7 and Table 3 summarise the information for amplitude, maximum and minimum, and first arrival for all thirteen cases.For the explicit model and the fine LEM model, the amplitude increases up to five times in an almost linear process as the stiffness of the fracture increases.The waveforms for the thick LEM end in an amplitude which is almost double that of the other two models.

Discussion
Mapping and understanding the creation and evolution of the EDZ is a topic of high importance when considering the safety of a future GDF.The P-wave velocity survey presented above is part of a series of experiments carried out in the future GDF of Finland in order to link the velocity changes with the man-made fractures (e.g., Refs.6, 10-13).However the objective of this work is to use methods and conclusions from previous studies 14,15,31,34 and assess them on field data.A source inversion method developed for laboratory studies 15,24 has been enhanced to create an input source for modelling a far more complex velocity survey, achieving a good result as shown by correlation coefficients ranging between 0.86 to 0.97.The flexibility of the LEM model in terms of thickness 15 led to two subcategories of the LEM approach one with thick and one with very fine LEM layers for each fracture set.The need to run large models in the most efficient way led to two approximate explicit fracture models, the pixelised PFM and the EDFM models.Modelling such a survey was a challenging process for many reasons, including the quality of the waveform data Table 3 Values of normalised peak-to-peak amplitude and values of normalised first arrival of for all cases used on optimising fracture stiffness for fracture F1 model waveforms.The values are normalised to the survey data.The peak-to-peak amplitude was calculated for the first part of the waveform with a time window of 0.4 ms (between time 0 to 0.4 ms of the recorded data).due to coupling, 13 the complex geometry of the fracture network, and the limited information on the stiffness of the fractures.The initial results using EM, two LEM (thick and fine) and two explicit (PFM and EDFM) models for representing fractures are in agreement with the conclusions from previous laboratory and theoretical studies. 15The full EM model did not work well in terms of arrivals, amplitudes, and full waveforms.The waveforms from the EM models are delayed 1.26 to 1.65 times compared to the survey waveforms (as a reference Pto S-wave velocity ratio is 1.72, Table 1).The reason is that the EM model has the same (slower) velocity for all frequencies (high as well as low frequencies) whereas the effect of fractures is to slow only low frequencies and to attenuate high frequencies. 18The amplitude of the EM waveforms is also 4 to 25 times higher than the recorded waveforms which is consistent with previous findings that the effective medium model cannot correctly predict the amplitude of transmitted waves, 15 (e.g., Refs.26, 35, 36).The frequency content is too high for the size of the fractures and the EM model is not expected to match the survey data, 15 but it is useful as a reference point for comparison of the two LEM cases.The EDFM model also fails to correctly predict the arrival time and the amplitude and the model waveforms are more attenuated than the recordings and delayed.The dip angles of the majority of the fractures are between 4 • to 18 • and these values are below the lower limit of 30 • above which the EDFM model produces more comparable waveforms. 15The LEM with thick layer works well for some cases on the initial models but it is not able to predict parts of the waveforms related to the reflections caused by surrounding fractures.Likewise, the pixelised explicit model (PFM) successfully matches the arrival time and amplitude in most cases in the initial models but it also underestimates the reflections requiring the stiffness of the EDZ fractures to be lowered.The fine LEM model is successful for the arrival time but not amplitude in the initial models.The manual iterative optimisation process for fracture stiffness has worked by increasing the amplitude of the waveforms and producing results which are more similar to the survey data for all three models.For the PFM and the fine LEM model, case 5 where the stiffness is   = 1670 GPa is the best case for uniform stiffness to match the survey data.The equivalent value for the thick LEM is   = 667 GPa for case 2. In comparing cross-hole tomography data against modelling 13 concluded that it was easy to match the first two pulse cycles but the latter part becomes more unpredictable.This statement is consistent with the current models.Although results have been improved, they never reach a full match with the survey data.One possible explanation is that the fracture stiffness used is uniform throughout the whole fracture.As the uniform stiffness increases the amplitude of the waveform crossing the fracture increases in a linear relationship.In nature, such a fracture does not exist.A stress-dependent stiffness will be very computationally intensive for such a large model, but an alternative solution would be to further optimise the fracture stiffness by creating manually zones with higher and lower stiffness in order to mimic the stiffness of a natural fracture and then automating the optimisation.

Conclusions
We have developed a methodology for modelling a seismic velocity survey which took place in the GDF ONKALO in Olkiluoto, Finland as part of a study for mapping the fractures in the EDZ.The purpose of this work was to test the performance of different fracture models in a real, complex engineering problem.The fractures of the studied area had been mapped and these were implemented within the models.Since the geometry of the mapped fractures is complex, their implementation within an orthogonal staggered grid needs to be approximated.The study compared five different methods of representing the fracturing.Two approaches were used for explicit representation of fractures based on Parastatidis 15 -the Pixelised Fracture Model (PFM), and the Equivalent Discrete Fracture Medium (EDFM).The fracturing was also represented by an Effective Medium (EM) and a Localised Effective Medium (LEM) 15 with two versions for the LEM model -the first using a thick LEM layer and the second using fine LEM layers.A wave source was derived for the models via an inversion process using waveforms from the deepest raypaths where there was no presence of fractures.A preliminary fracture stiffness was used based on the size of each fracture following Worthington and Lubbe 32 and this was then optimised based on waveform comparisons.Due to the vast number of waveforms (552 waveforms per survey), comparisons were made for selected raypaths based on different levels of fracture influence.In the first phase of modelling, neither the EM model nor the EDFM model compared well with the data.The EM model had much higher amplitudes and late arrivals while the EDFM had over attenuated amplitudes and late arrivals.The explicit PFM model and the two LEM models compared more favourably with the data and hence a second phase which optimised the fracture stiffness for the F1 fracture was conducted for these models.The main conclusions of this investigation are: • The two LEM models and the explicit PFM model worked very well for some cases in terms of reproducing velocities, amplitude and frequencies.However, further optimisation of the fracture stiffness will be needed to realistically reproduce the waveforms.• The EM model cannot give a reasonable solution to a problem with such a complex fracture network as the EDZ.Averaging all fracturing into a single medium removes all heterogeneity and in particular results in unrealistically high amplitudes.• The EDFM model needs an improved approach to balance the size of the equivalent discrete fractures with the fracture density and the fracture stiffness in order to achieve more realistic results.• Using the LEM with thick layers for each fracture can give a quick and easy-to-build model (no need to write model files We have selected and presented 3 cases out of the 13: Case 2 (  twenty times higher than initial), Case 4 (  forty times higher than initial), and Case 5 (  fifty times higher than initial).
with complex objects) with reasonably good results.The resulting waveforms are closer to those of the survey than those produced by the EM model.The results however are not as good as the explicit PFM model or the fine LEM model.• Optimising fracture stiffness for a single fracture based on a single source-receiver path improved results for the three models, increasing the amplitudes such that the first two pulses of survey and modelled data are comparable.• The thick layer LEM model was more sensitive to changes in fracture stiffness and it required lower values of fracture stiffness than the explicit PFM model and the fine LEM model to obtain the correct survey amplitudes.In general, it produced higher amplitudes and this is clearly related to the problems with the homogeneous EM model which produces unrealistically high amplitudes.• As found in previous work, the explicit model and the LEM model with fine layering perform similarly and they produce the same optimal values for fracture stiffness.
Overall, none of the models are able to match the survey results perfectly.This is not simply due to geometrical approximation of the fractures but is a limitation of assuming a uniform fracture stiffness for the fracture in the optimisation process.Previous studies have shown that for fractures that are large relative to wavelength the variation of fracture stiffness across the fracture significantly influences the wave behaviour and needs to be encapsulated in the model. 14,18For this reason further optimisation needs to be made to the fracture stiffness either by creating zones with different stiffness in the fracture or by applying stress-dependent fracture stiffness.For the latter, further code development will be required.

CRediT authorship contribution statement
Emmanouil Parastatidis: Conceptualization of this study, Methodology, Software.
where h is the thickness of the LEM layer.A way to convert crack density  to cracks per unit length 1  and vice versa is then Eq. ( 8) presents the stiffness matrix calculated from five independent material constants linking cracks per unit length 1/L, fracture stiffness (  ,   ), and 2 elastic constants for an isotropic background material with derived parameters shown in Eq. ( 9), Coates and Schoenberg.
where  is Lamé's first parameter,  is Poisson's ratio and  is the shear modulus (Lamé's second parameter).

B.2. Examples of manual iterative optimisation for more raypaths
See Figs.11 and 12.

Appendix C
See Table 4.We have selected and presented 3 cases out of the 13: Case 2 (  twenty times higher than initial), Case 4 (  forty times higher than initial), and Case 5 (  fifty times higher than initial).We have selected and presented 3 cases out of the 13: Case 2 (  twenty times higher than initial), Case 4 (  forty times higher than initial), and Case 5 (  fifty times higher than initial).

Fig. 1 .
Fig. 1.Position of testing tunnel niche ONK-TKU3620 in the ONKALO Underground facility where the EDZ studies where conducted at 345 m depth.16The image shows the extent of ONKALO that was available during the execution of the EDZ experiment, before being expanded into GDF.

Fig. 2 .
Fig. 2. Schematic view of the EDZ and how is divided based on the origin of the fractures in stress induced (EDZ  ) and construction induced (EDZ  ) and EdZ using the drilling and blasting method (D&B) used in Olkiluoto as presented in Ref. 16 and modified after.2

2
Fig. 2. Schematic view of the EDZ and how is divided based on the origin of the fractures in stress induced (EDZ  ) and construction induced (EDZ  ) and EdZ using the drilling and blasting method (D&B) used in Olkiluoto as presented in Ref. 16 and modified after.2

Fig. 3 .
Fig. 3. (a) The 3D lithological model of the whole testing tunnel with the contact area between veined gneiss (VGN) and pegmatitic granite (PGR).The zoomed area is the EDZ study highlighting the modelled area and E4 block (modified from Refs. 7, 19) (b) 2D projection of the testing tunnel ONK-TKU-3620 niche that the velocity survey area took place and the modelled area E4, 7,8 the red marked area (blocks E4-E7) is the sawed blocks, the velocity survey took part in blocks E3-E8 (c) Fracture model based on the wire sawing results (blocks E4-E7) and extrapolated to cover the whole area based on available drillhole investigations (modified from Ref. 20), the red fracture is stress-induced (F1) and boreholes 34-43 are the modelled area.(d) One 2D tomography plane consists of 4 surveys with 24 shots.The red circles are the position of the receiver/transducer, and the blue lines are the raypaths created when transducing from each of these positions. 8(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

E
.Parastatidis et al.

Fig. 4 .
Fig. 4. (a) Approximate 2D cross-section of the fracture model for borehole frame 34-43 for E4 block to be modelled.(b) Cross-section of the PFM, (c) the EDFM model, (d) the EM fracture model, (e) the thick layer LEM fracture model and (f) the fine layer LEM fracture model for the borehole frame 34-43, fracture F1 (blue) number N 7 in Appendix C Table 4 and EDZ fractures (red) number N 9, 11, 12 and 15. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Waveforms for the frame between boreholes 34-43 for the 5 models.The shot position is on borehole 43 and recording on borehole 34.5(a), Sketch of the raypath, for the shot depth at 0.6 m and recording at 0.55 m, 5(b) PFM model, 5(b) EMDF model, 5(d) EM model, 5(e) thick LEM and 5(f) fine LEM.

Fig. 6 .
Fig. 6.Values used for optimising Normal   and Shear   fracture stiffness for fracture F1 versus model case.

Fig. 7 .
Fig. 7. Peak-to-peak amplitude versus the normal 7(a) and shear 7(b) stiffness values for all 13 cases of stiffness optimisation (Fig. 6), the amplitude values are normalised, where 1 is the maximum peak-to-peak amplitude (for the first 0.4 ms) of the survey data.

Fig. 9 .
Fig. 9. Waveforms for the frame between boreholes 34-43 for the 5 models.The shot position is on borehole 43 and recording on borehole 34.On the side is the sketch of the raypath of the wave 9(a), shot depth at 0.6 m and recording at 0.20 m, 9(b) PFM, 9(b) EMDF models, 9(d) EM model, 9(e) thick LEM and 9(f) fine LEM.

Fig. 10 .
Fig. 10.Waveforms for the frame between boreholes 34-43 for the 5 models.The shot position is on borehole 43 and recording on borehole 34.On the side is the sketch of the raypath of the wave 10(a), shot depth at 0.6 m and recording at 0.60 m, 10(b) PFM, 10(b) EMDF 10(d) EM model, 10(e) thick LEM and 10(f) fine LEM.