Investigating the use of 3-D full-waveform inversion to characterize the host rock at a geological disposal site

The U.K. government has a policy to dispose of higher activity radioactive waste in a geological disposal facility (GDF), which will have multiple protective barriers to keep the waste isolated and to ensure no harmful quantities of radioactivity are able to reach the surface. Currently no speciﬁc GDF site in the United Kingdom has been chosen but, once it has, the site is likely to be investigated using seismic methods. In this study, we explore whether 3-D full-waveform inversion (FWI) of seismic data can be used to map changes in physical properties caused by the construction of the site, speciﬁcally tunnel-induced fracturing. We have built a synthetic model for a GDF located in granite at 1000 m depth below the surface, since granite is one of the candidate host rocks due to its high strength and low permeability and the GDF could be located at such a depth. We use an effective medium model to predict changes in P -wave velocity associated with tunnel-induced fracturing, within the spatial limits of an excavated disturbed zone (EdZ), modelled here as an increase in fracture density around the tunnel. We then generate synthetic seismic data using a number of different experimental geometries to investigate how they affect the performance of FWI in recovering subsurface P -wave velocity structure. We ﬁnd that the location and velocity of the EdZ are recovered well, especially when data recorded on tunnel receivers are included in the inversion. Our ﬁndings show that 3-D FWI could be a useful tool for characterizing the subsurface and changes in fracture properties caused during construction, and make a suite of suggestions on how to proceed once a potential GDF site has been identiﬁed and the geological setting is known.


I N T RO D U C T I O N
Over the last 70 yr, a large legacy of radioactive waste has accumulated in the United Kingdom. A significant amount of higher activity waste (HAW) has been accrued (NDA 2015) and now needs to be securely isolated from the surface biosphere. To ensure the HAW is safely contained over geological timescales (>100 000 yr), it will be disposed in a deep geological disposal facility (GDF), to be built at a depth of 200-1000 m below surface. The GDF will consist of multiple components including engineered, chemical and geological barriers and will have a facility footprint that could be about 10 km 2 (though the associated surface facilities will be much smaller in extent). At the moment, no site has been selected, but several potential host rocks have been identified including 'soft' rocks (e.g. clays and mudstones), 'hard' rocks (e.g. granite) and halite.
Granite is considered a potentially suitable host rock for radioactive waste disposal because it has high bulk strength and low ground water permeability, and several countries have already undertaken geophysical investigations to site GDFs in granitic rocks, for example Finland (Cosma & Heikkinen 1996;Saksa et al. 2007;Schmelzbach et al. 2007;Cosma et al. 2008) and Sweden (Juhlin et al. 2002;Bergman et al. 2006;Juhlin et al. 2010). Though fluid flow and radionuclide transport in granites are reasonably well understood, it is challenging and necessary to assess how rock permeability may change during and after the construction of tunnels (Jaeger et al. 2007). Typically, tunnelling in granite produces two distinct regions under stress: the region closest to the tunnel, the excavated damage zone (EDZ), which is subject to irreversible damage; and the next radial region, the excavated disturbed zone (EdZ), 1 where the changes are elastic and recoverable (Tsang et al. 2005). The EDZ (which may be 2-3 m thick) is the most fractured, and its characteristic properties include stress anisotropy, changes in fracture density and orientation, and enhanced or reduced fluid flow. Additionally, it is difficult to reduce fluid flow through this zone by sealing (grouting) after it is damaged due to the back pressure from the rock mass (Tsang et al. 2005).
In contrast to the EDZ, the EdZ is defined as a region where only reversible elastic deformation occurs (Tsang et al. 2005). Even though no new fractures are formed in the EdZ, the changes in stress field can temporarily alter the properties of the existing fracture network, such as opening existing fractures, for an unknown length of time. The EdZ may extend to large distances (10 s of metres) from the tunnel (Perras & Diederichs 2016) but it has not been possible to accurately define the outer limits of the zone (Tsang et al. 2005). It can be difficult to detect and monitor the presence of an EdZ in situ; for example, at theÄspö Hard Rock Laboratory (Sweden) where several tunnel-based surveys failed to either detect the EdZ or evaluate its fracture properties (Siren et al. 2015). Consequently, improving technology to detect the EdZ and monitor changes in its hydromechanical and geochemical processes, is essential for the long-term safety of a GDF (Tsang et al. 2005).
To characterize the EdZ and host rock in situ, seismic imaging can be used (e.g. Cosma & Heikkinen 1996;Hildyard & Young 2002;Pettitt et al. 2006;Schmelzbach et al. 2007;Juhlin et al. 2008;Marelli et al. 2010;Zhang & Juhlin 2014;Reyes-Montes et al. 2015) since seismic properties are influenced by rock fracture patterns and properties. In particular, surface seismic reflection surveys are useful in characterizing host rock bulk properties and large-scale structures such as faults that may affect the geological barrier integrity of a GDF (Cosma & Heikkinen 1996). Additionally, seismic reflection surveys in facility tunnels have produced detailed images of fracture networks (Cosma et al. 2013;Brodic et al. 2017). Furthermore, 2-D full-waveform inversion (FWI; Zhang & Juhlin 2014) and microseismic/acoustic emissions (AE) methods (King et al. 2011;Saari & Malm 2013;Goodfellow & Young 2014;Reyes-Montes et al. 2015) have been successful in mapping fractures in proposed GDF host rock, and tunnel surface waves generated during tunnelling have been used to detect geological structures ahead of the tunnel (Jetschny et al. 2011).
In this study, we investigate the use of seismic data to characterise a potential GDF site with a hypothetical stress induced EdZ associated with the construction of tunnels in a granitic host rock. Building on the successful application of 2-D FWI by Zhang & Juhlin (2014), we explore whether commercially used 3-D FWI codes with additional capabilities (Warner et al. 2013;Debens 2015;Warner & Guasch 2016;Agudo et al. 2018a,b,c) can improve the detection of the EdZ in situ around a model GDF. Additionally, we start with a conventional surface array, and investigate what source-receiver offsets are required to recover the subsurface structure well. We then add receivers within the tunnel to see whether this enhances the resolution of the EdZ. Finally, we consider the effect of tunnel infrastructure on the performance of FWI and our innovative survey designs in this particular granitic GDF environment.
Using our workflow, we first build a numerical model of a hypothetical granitic GDF environment and assign representative P-wave velocities. In the model the EdZ is placed within the host rock and is characterized by reduced P-wave velocity, caused by increased fracturing. Additionally, we develop and test more complex models that contain tunnel infrastructure. Next, we generate seismic data for each of our velocity models using our two different survey designs: a conventional surface-survey; and a combined survey with sources and receivers at surface, and receivers within the tunnel.
We apply 3-D FWI to recover a model of velocity across the EdZ using a starting model in which P-wave velocity values are those of undamaged rock, that is without the tunnel-induced EdZ. Then, we assess the recovery of the inverted EdZ target for both survey designs and evaluate the effect of survey size and tunnels on the overall inversion result. In general, we find that 3-D FWI is successful in resolving the EdZ in our selected granitic host rock and, in particular, we discover that the combined survey is important for good recovery of the EdZ for seismic surveys with reduced shotreceiver offsets and for models that include tunnel infrastructure. Finally, we suggest some further tests that could be performed once potential locations for the GDF site have been identified.

Geological setting-granitic host rock with sedimentary overburden
We select a granitic host rock for our hypothetical GDF. The model consists of an 800-m-thick sedimentary overburden, and a 400m layer of fractured granite above unfractured granite bedrock, as this lithological combination is found in the United Kingdom (e.g. Towler et al. 2008). The GDF is located in the middle of the fractured granite at a depth of 1000 m. In addition, the model contains an anomalous zone around the GDF. During construction, tunnelling could locally increase fracturing and/or open pre-existing fractures, resulting in EDZ and EdZ in the tunnel walls. For our initial inversions we keep the velocity model simple, and define a single large combined EdZ (individual shafts and tunnels are omitted) characterized as a low-velocity zone caused by an increase in fracture density due to tunnelling (see Section 2.2). The EDZ is not explicitly modelled since it is too small (2-3 m) to be detected using the modelling and seismic survey parameters chosen in this study (see Sections 3.1.2 and 3.2). The dimensions of the modelled EdZ are: 500 m x 500 m horizontally, well within the expected footprint of the underground facility and 150 m vertically, to account for a potential facility design with tunnels at multiple depths. 2 We note that, in order to model smaller more realistic target features, we would have to decrease the grid spacing in the velocity model, use a denser array of shots and receivers, and compute the wavefield with a reduced time step (see Section 3), which all significantly increase the computational effort. For our preliminary tests we are principally interested in exploring whether FWI can resolve an EdZ, so have used the same EdZ anomaly size in all the models shown here (see Section 4).

Medium properties
We assign geophysical properties for the geological model described in Section 2.1 that are consistent with typical values. Most notably, P-wave velocity is larger in the granitic basement than in the sediments. Additionally, we calculate gradual increases in velocity with depth within each unit to reflect compaction effects (Barton 2006). To determine changes in P-wave velocity in the host rock due to increases in fracture density, S-wave velocity and density are also required for this layer (Table 1). Table 1. Key parameters used to define the three lithological units and EdZ in the true velocity model. P-wave velocity (Vp) increases within each unit linearly with gradient of 1 or 0.46 m 2 s -1 for the granitic units (Barton 2006). Reference isotropic parameters for granite host rock at 1000 m: Vp = 5850 m s -1 ; Vs = 3400 m s -1 ; ρ = 2850 g cm -3 .
We calculate the influence of vertical fractures on the P-wave velocity using Equivalent Medium Representation (otherwise known as Effective Medium theory) of Horizontal Transverse Isotropic (HTI) media. More specifically, we model the vertical fractures as small isolated circle cracks in a planar distribution to find the effective P-and S-wave velocities and density (Liu et al. 2000). To implement this model, we choose to define the fracture density in the fractured host rock as 0.1 and double the fracture density to 0.2 in the EdZ, and for all the units, we calculate the 6 × 6 elastic stiffness matrix for HTI media. As expected, elastic constants C 11 and C 33 , seismic wave velocities and rock density all decrease when fracture density increases (Fig. 1). Since the change in velocity follows an exponential decrease, the largest changes in velocity are observed for fracture densities of <0.2. It should be noted that the fracture density of the host rock and EdZ could vary at these depths and fractures could be dry, wet or closed. Changes in fracture density will likely yield similar changes in seismic velocity as shown in our models (Fig. 1). Additionally, substituting properties for wet fractures should produce velocity anomalies of a similar order of magnitude to dry fractures (Liu et al. 2000). As such, the workflow developed here is appropriate for such changes in velocity.

Imaging challenges
Given the geological setting selected for our tests, with a highvelocity granitic rock overlain by a comparatively lower velocity sedimentary layer, there are imaging challenges in characterizing the host rock and EdZ. The most exploited seismic phase in nearsurface exploration, the reflected P wave, has a small critical angle in this setting. Most of the energy of the wave is reflected rather than transmitted, and thus is more suited to imaging changes in reflective coefficients than velocity structure. In contrast, seismic refraction (transmission) waves are more sensitive to medium-tolong wavelength velocities (Pratt et al. 1996;Sirgue 2006;Vireux & Operto 2009) and, in our case, are useful in revealing the granitic host rock velocity structure. For our study, we include relatively long-period transmission waves and source-receiver offsets that are large enough to ensure the seismic wavefield passes through the host rock and EdZ.

FWI overview
FWI is a computational technique for generating high-resolution, high-fidelity models of physical properties in the subsurface. It is a local, iterative inversion scheme that successively improves a starting model. It uses the two-way wave equation to predict seismic data from the starting velocity model, and updates this model in a way that minimizes the difference between the predicted and observed data (Warner et al. 2013).
The use of FWI has expanded rapidly in the last 10-20 yr, and many industry and academic groups have developed their own software. The most important advance for the petroleum sector was the move from a 2-D to 3-D scheme and, only then, was FWI considered to be of commercial use (Sirgue et al. 2010;Warner et al. 2013;Operto et al. 2015). 2-D FWI can recover accurate velocity models, but 3-D FWI leads to improved recovery (Agudo et al. 2018b), even for seismic profiles that are close to 2-D (Kalinicheva et al. 2017). The next most significant advancement was the addition of anisotropy, which means that the kinematics (traveltimes) of the wavefield can be correctly predicted. 3-D acoustic, anisotropic FWI has now been widely adopted by the petroleum sector, and has been demonstrated to be successful using advanced quality assurance procedures, including noticeable improvements in 3-D pre-stack depth migration images, improved flattening of common image gathers and better matches with borehole data (Prieux et al. 2011;Kapoor, et al. 2013;Selwood et al. 2013;Warner et al. 2013). The geological disposal site will be complex in three dimensions; most near-surface rocks are anisotropic and any induced fracturing will produce additional anisotropy, hence, the ultimate use of a 3-D FWI code with anisotropy is warranted.
The code utilized here can solve for tilted transverse isotropy (TTI) anisotropy which, as described above, will be an important capability in any application of FWI to the field dataset acquired across the GDF. The code also has the capability of alternating between FWI (a local inversion) and a global inversion for anisotropy and attenuation (Debens 2015;da Silva et al. 2017). If anisotropy is not accounted for, FWI velocity models are stretched and the depths are inaccurate, as seen in the recent drilling of the Chicxulub impact crater, in which faster subhorizontal FWI-determined velocities in the sedimentary overburden led to an overestimation of depth to the crater (Christeson et al. 2018). We have not included anisotropy in the tests shown here as this approximately doubles the computational effort (Warner et al. 2013).
There is also an option to model and/or invert for the elastic wavefield, but it is also computationally very expensive, and is rarely required in marine data sets. We have encountered one single case where the acoustic code failed due to not adequately accounting for the elastic properties of the wavefield (Agudo 2018;Stronge 2018). After the (elastic) field data were converted to acoustic data using a Wiener filter matching scheme (Agudo et al. 2018a), however, an acoustic inversion of the matched data was successful (Agudo 2018).
A common problem with performing FWI is cycle skipping, which leads to a recovered velocity model that is located in a local rather than global minimum (Pratt 1999). Cycle skipping occurs if the starting model is unable to predict the majority of data to within half a cycle of the field data (Sirgue et al. 2010) at the lowest inversion frequency. This means that, typically, significant effort has to be expended on obtaining a starting velocity model that is close to the true model. To address this issue a second code was developed with a different objective function, adaptive waveform inversion (AWI), which is less sensitive to cycle skipping (Warner & Guasch 2016;Guasch et al. 2018;Yao 2018). We use AWI for initial inversions when the starting model is poor and then move back to FWI once the recovered velocity model has improved sufficiently such that it is no longer cycle-skipped. Effectively this means it is possible to start with a poor starting model and still get to the correct answer, albeit after a larger number of iterations.
For this study, where we invert for synthetic data, we use the 3-D acoustic inversion scheme for computational efficiency. The algorithm proceeds as follows: (1) Calculate the direction of the local gradient (i) Using the starting model and a known source, calculate the forward wavefield everywhere in the model including at the receivers.
(ii) At the receivers, subtract the observed from the calculated data to obtain the residual data.
(iii) Treating the receivers as virtual sources, back-propagate the residual data into the model, to generate the residual wavefield.
(iv) Scale the residual wavefield by the local slowness, and differentiate it twice in time.
(v) At every point in the model, cross-correlate the forward and scaled residual wavefields, and take the zero lag in time to generate the gradient for one source.
(vi) Do this for every source, and stack together the results to make the global gradient.
(2) Find the step length (i) Take a small step directly downhill, and recalculate the residual data.
(ii) Assume a linear relationship between changes in the model and changes in the residual data.
(iii) Use the resulting straight line to decide how far to step downhill to reduce the residual data to zero.
(iv) Step downhill by the required amount, and update the model. This procedure is repeated until changes to the model become minimal. The computational effort required for FWI is large, but the resulting spatial resolution is much better than can be obtained by methods that seek to match traveltimes, for example first-arrival traveltime or reflection tomography.

Implementation of FWI
Using the full seismic waveform for the velocity inversion gives us the potential to detect subtle changes in medium properties and structure in the host rock. To apply FWI in this setting we simulate surface and combined seismic surveys, and use maximum offsets of around three to five times the depth of the GDF to ensure that the transmitted wavefield passes through the target. Following the workflow identified in Section 3.1.1, synthetic data are generated for each velocity model and used as the 'observed' input data for FWI. We assume that the background velocity is known reasonably well, and use a starting model that is identical to the true model, except without the EdZ. This is justified for this hypothetical case but for real data we would ultilize AWI (as described in Section 3.1.1), which is less sensitive to cycle skipping and allows us to obtain a good starting model for FWI from a relatively poor starting model (Warner & Guasch 2016;Guasch et al. 2018). For computational efficiency, the acoustic wave equation is used and thus only a P-wave velocity field is required (Table 1). We note that, the vast majority of industry applications of 3-D FWI use an acoustic approximation and are successful (Vireux & Operto 2009;Bansal et al. 2013;Kapoor et al. 2013;Selwood et al. 2013;Warner et al. 2013;Operto et al. 2015). We adopt the multiscale approach that is widely used in FWI applications, by gradually inverting data with an increasingly higher frequency content. We start our inversions by inputting data up to 8 Hz, as this is the lowest frequency able to detect the spatial features in the true velocity model.
For numerical modelling, we need to define a suitable grid structure and time sampling that ensures stability and limited dispersion. We selected a grid spacing (dx) of 12.5 m since this allows us to include tunnels in our synthetic velocity model. The maximum frequency (f max ) we can model with this grid spacing is: where P-wave velocity (Vp min ) is 3700 m s -1 and number of points per wavelength (n) is 4. Therefore, the maximum stable frequency is 74 Hz. We choose a sufficiently small time step such that the wavefield travels no more than half a grid cell in a single time step. As the highest rock velocity is 6150 m s -1 , we use a time step of 1 ms. It should be noted, however, that in the models with tunnel infrastructure, the minimum P-wave velocity is 342 m s -1 (in the tunnels) and the maximum inversion frequency should thus be set to 6.8 Hz to avoid dispersion and modelling artefacts. We are able to use higher frequencies when inverting synthetic data, but we note that we would have to use a smaller grid spacing if we wished to invert for frequencies of up to 74 Hz when applying FWI to field data.
Our early tests indicate that we are able to recover the velocity anomalies reasonably well using 32 iterations across six frequency bands, with four iterations for each band and maximum frequencies of 8, 12 and 17 Hz, and eight iterations for bands with maximum frequencies of 24, 33 and 43 Hz. We use the exact same number of iterations and inversion frequencies for all the results presented here to evaluate the performance of FWI for different geological models and survey designs. Even though a maximum frequency of 74 Hz can be used, we find this is not necessary to recover the structures in the model. Boundary conditions are applied to each of the six model boundaries. For the top boundary, we assign a free surface condition, allowing the energy to reflect back into the model space.
At the sides and bottom of the model, we use absorbing boundary conditions. We assess the quality of the inversion through analysing the difference between the inverted and true velocity field (Section 4). Additionally, we check that the global functional decreases with increasing iterations within a frequency band.

Sources and receivers
A key interest of this study is to explore whether combined surveys utilizing surface and tunnel receivers, improve the inversion result in comparison to conventional surface surveys. For the surface surveys, we use 3861 surface receivers with a spacing of 50 m buried at 12.5 m depth (Fig. 3). For the combined surveys, we use the same surface receiver geometry and include either: For all surveys, we use 250 surface sources at a depth of 12.5 m (Fig. 3). All receivers and sources are located on a regular grid (Fig. 3).

Coverage found using ray tracing
Prior to generating synthetic data, we use raytracing to assess the coverage of turning and reflected waves within the target EdZ. This fractured zone is anisotropic, therefore we use the fully anisotropic (up to triclinic symmetry) seismic ray tracer ATRAK (Guest & Kendall 1993) and implement the velocity and density fields as listed in Table 1. For many take-off angles (and therefore sourcereceiver offsets) there is substantial ray coverage within the EdZ. For example, for take-off angles 30-60 • , turning rays travel through the lower section of layer 1 (sedimentary overburden) and through the EdZ in the centre of the model (Fig. 4). Additionally, for these angles, rays that reflect at the interface of layers 2 and 3 (i.e. the base of the fractured granite host rock) travel through the EdZ allowing increased and diverse sampling through the target region. Further testing using raytracing confirms that velocity gradients are essential in producing turning waves that travel through the EdZ.

R E S U LT S
Using the geometry defined in Section 3.2 we analyse the effectiveness of combined surveys and 3-D FWI in resolving velocity structures for three different velocity models. Model 1 is a simple three-layer model and contains the three lithological units and the 500 m x 500 m x 150 m EdZ, the structural target for the study. Models 2 and 3 are based on Model 1 but also include some tunnel infrastructure. Model 2 has a simple tunnel layout with one disposal tunnel at 1000-m depth and one vertical access tunnel connecting disposal tunnel and the surface. Model 3 has a more realistic tunnel arrangement with six disposal tunnels at 1000-m depth, as well as the vertical access tunnel. The tunnels are 25-m wide and have a P-wave velocity of air (342 m s -1 ). In the following sections, we show the results from the simple model (Model 1) and increase the model complexity through adding tunnels (Models 2 and 3) so we can evaluate whether a combined survey is important in resolving the disturbed zone.

Velocity models
The simple three-layer velocity model and EdZ have properties as listed in Table 1 and is shown in Fig. 2. The starting model used to initiate the inversion contains the same lithological layers as the true velocity model (Fig. 2) but does not include the low velocity disturbed zone.

Results
The resultant velocity models (with starting model subtracted) show that both the surface and combined surveys resolve the EdZ well (Fig. 5). We next investigate how critical the long-offset data are for the inversion, in consideration of a scenario whereby the extent of the survey area is limited. We evaluate the performance of surveys with a reduction in both survey area and maximum offset. We quantitatively compare the inverted velocity fields by computing a model misfit, defined as the total RMS misfit between the inverted and true model profiles. We look at trends in the RMS misfit for three locations: outside, at the edge and in the centre of the EdZ (Fig. 6). For both survey types, the total RMS misfit increases when decreasing the survey area and maximum offset, as expected. Additionally, for all maximum offsets for all three locations, the combined survey has lower RMS misfits than the surface survey. In the centre of the EdZ, the largest improvements in inversion result are observed since the combined survey has 50 per cent lower misfit than the surface survey ( Fig. 6). At the edge of the EdZ, the difference in RMS misfit is reduced when increasing the maximum offset (survey area size) such that no significant difference is observed using 5 km offsets. Furthermore, we observe that including the tunnel receivers in the inversion reduces the misfit when the survey area is restricted.

Velocity models, tunnel infrastructure and properties
As shown in Section 4.1, the EdZ is detectable when using either the surface and combined surveys and the maximum offset range. To further test these survey designs we consider including complexity in the velocity model by adding GDF tunnel infrastructure (Figs 7a  and b). We choose the tunnels to be 25-m-wide open cavities, 3 with P-wave velocity equivalent to wave speed in air (342 m s -1 ). During the inversion, the tunnel velocity is fixed and not updated, since the location and dimensions of the tunnel are assumed to be known.
The tunnel system design is reasonably basic with two orthogonal tunnels: a horizontal disposal gallery (Fig. 7b) and a vertical access shaft (Fig. 7a). The EdZ remains the same size, 500 m x 500 m, for consistency (Fig. 7b), though we appreciate that the EdZ for a single tunnel would be smaller (10 s of metres). Likewise, the 36 tunnel receivers will have the same geometry as in Section 4.1, despite any practical acquisition restrictions.

Results
In Fig. 8 we subtract the starting velocity field from the true velocity field; the inverted velocity using surface survey; and the inverted velocity using the combined survey to allow focus on the EdZ low velocity zone. For most depths, there are minor differences between the inverted fields. The largest variations in recovered velocity are seen for depths below the tunnel receivers (>1050 m), especially at 1062.5-m depth (Fig. 8). Although both surveys detect a velocity anomaly associated with EdZ, for the surface survey the shape of the low velocity region could be interpreted as two separate anomalies, due to the disposal tunnel creating a poorly resolved region (Fig. 8b). By incorporating the tunnel receivers into the combined survey, the shape of the EdZ is resolved more completely (Fig. 8c) and more closely resembles the true anomaly (Fig. 8a).

Velocity models, tunnel infrastructure and properties
As shown in Section 4.2, the combined survey resolves the EdZ more completely than using surface receivers only. We advance the tests by increasing the complexity of the GDF tunnel infrastructure by implementing five parallel horizontal disposal galleries at 1000m depth (connected by a horizontal access tunnel) in additional to the vertical access shaft. The tunnels remain as 25-m-wide open cavities with P-wave velocity equivalent to wave speed in air (342 m s -1 ). The extent of the EdZ is not only consistent with the previous model but is now an appropriate size for a combined EdZ in complex tunnel system. As the Model 3 contains more tunnels and potentially more poorly resolved regions, we compare two configurations of tunnel receivers. The first tunnel survey is the same as the survey used in Sections 4.1 and 4.2 with 36 tunnel receivers at 1050 m depth and separated by 100 m. The second tunnel survey analysed has 50 receivers in total at two depths different 950 and 1050 m, essentially 50 m above and below the centre of the disposal tunnels, and are also separated by 100 m.

Results
Similar to Section 4.3.1, we complete the inversions for both surveys using a true velocity model with complex tunnel infrastructure. Again, to focus on how well we resolve the EdZ low velocity zone, the starting velocity field is subtracted from the true velocity model; the inverted velocity field using the surface survey; and the inverted velocity fields using the two combined survey layouts. The largest Downloaded from https://academic.oup.com/gji/article-abstract/215/3/2035/5101440 by University of Leeds user on 05 December 2018 variation in wavefield sampling occurs above and below the tunnels (due to the complex tunnel footprint), and thus the results are shown in cross-section to highlight the key differences between the survey designs ( Fig. 9).
For most depths (0-900 and 1150-2000 m) there are minor differences between true and the inverted fields (Fig. 9). However, there are distinct differences within the resolved EdZ, caused by receiver geometries. The surface survey reasonably resolves the EdZ above the tunnels but does not completely recover the velocity structure beneath the tunnels (Fig. 9b). The first combined survey (with receivers below the tunnels) improves the recovered velocity field below the tunnels (Fig. 9c). The second combined survey (with receivers above and below the tunnels) has the best recovery of the shape of the velocity anomaly above and below the tunnels (Fig. 9d). Additionally, the size of the velocity anomaly is recovered well at the edges of the tunnels at 950-1050 m deep, but velocity anomalies above and below the tunnels are not fully recovered (Fig. 9d).

D I S C U S S I O N
The synthetic tests presented here suggest that 3-D FWI may be a useful tool in characterizing a potential GDF site and detecting changes in rock properties associated with tunnel construction. We note, however, that the modelling is quite simplistic and that there are additional challenges associated with inverting a real field data set. The GDF site used here is purely hypothetical and may, or may not, be a good analogue for the future site. Once potential sites are identified by the Radioactive Waste Management (RWM), they will almost certainly commission the acquisition of seismic data to characterize the site. With this in mind, we recommend that further synthetic tests be performed to help design a base and any future seismic surveys, and ascertain whether 3-D FWI will be able to characterize the site, before and after tunnel construction.
With regards to the design of the future seismic survey, in general it is preferential to have randomly spaced shots and receivers rather than positioning the array on a regular grid. Regular grids tend to lead to linear artefacts along the grid lines (Warner et al. 2013). In addition, as shown here in Section 4.1, tests should be performed to determine what shot-receiver offsets are required to obtain refracted arrivals that penetrate the proposed depth of the GDF facility, which should improve the performance of FWI. Although reflections can be included in FWI (e.g. Yao et al. 2018), wide-angle refractions are important for recovering the medium-to-long wavelength velocity structure (Pratt et al. 1996;Sirgue 2006;Virieux & Operto 2009). Furthermore, we anticipate that the receiver spacing should be set to be approximately equal to the smallest anomaly size that RWM wish to resolve (Morgan et al. 2016).
It is appreciated that the number of receivers used in a surface survey is dependent on the size of the survey area available and as such we have shown that the inversion procedure does reveal the EdZ even with a small number of receivers and reduced maximum offset (Section 4.1). In practice, challenges in survey area such as topography can be overcome with new wireless technology (e.g. Savazzi & Spagnolini 2008;Crice et al. 2015). Such advances are useful in tunnel surveys too, but for this study the number and distribution of tunnel receivers deployed is quite conservative (Section 3.2.1, Fig. 3). Nevertheless, we demonstrate that the more distributed receivers are across the tunnel networks, the better the recovery of the EdZ (Sections 4.2 and 4.3, Figs 8 and 9). Additionally, preliminary inversion tests show that the EdZ is, perhaps unsurprisingly, better recovered when using more tunnel receivers in total (Section 4.3, Fig. 9).
The synthetic FWI tests presented here could be re-run using a newly constructed GDF model that matches the geology of the selected site. Additional inversions are also recommended, in order to make the simulations more representative of a FWI application to a real field dataset. For example, in Morgan et al. (2013), random noise was added to the synthetic data, the inversions were started with less accurate starting velocity models, the synthetic data were generated with an elastic code, and windowed in time so that only the first-arriving refractions were allowed into the inversions. Data windowing is often applied to field data prior to input to FWI, with short-offset reflections and secondary arrivals being removed through muting (e.g. Warner et al. 2013). Performing more realistic inversions will provide confirmation as to whether it is possible to recover the quite small velocity anomalies induced by tunneling and, perhaps more importantly, what coverage (experimental geometry) is needed to do so.
In the modelling shown here, we have used a starting model that has accurate background velocities. In future tests, the starting velocity model could be obtained through a traveltime tomographic inversion of synthetic data acquired across the new GDF model. For many years, the success of FWI has been dependent on having a good starting model, which needs to be able to generate synthetic data that are not cycle-skipped with the observed data (Pratt 1999;Sirgue 2006;Warner et al. 2013). Two approaches that can mitigate problems with poor starting models when performing 3-D FWI are the use of: (1) phase plots to identify and remove cycle-skipped data (Shah et al. 2012;Warner et al. 2013;Morgan et al. 2016) and (2) AWI for the initial iterations until the synthetic data are not cycle-skipped (Warner & Guasch 2016). Other approaches that have been developed to address cycle-skipping include: Optimal transport (Metivier et al. 2016), dynamic warping (Ma & Hale 2013) and tomographic FWI (Biondi & Almomin 2012). These schemes mean that it is possible to start with a relatively poor starting models and still recover the true velocity model.
With regards to using an acoustic rather than elastic code. Acoustic 3-D FWI codes have been successfully applied to many marine data sets, and their use is now standard practice within the petroleum sector (Bansal et al. 2013;Kapoor et al. 2013;Selwood et al. 2013;Warner et al. 2013;Operto et al. 2015). Any future seismic data acquired across a land GDF facility are, however, likely to be more strongly affected by elastic effects. One scheme that could be used here is a transformation of the elastic (field) data to acoustic data, which has been shown to improve the recovery of the true velocity structure using an acoustic inversion, for both marine and land data (Agudo et al. 2018a). In addition, deriving an accurate source is more challenging for land seismic surveys (Rowse & Tinkle 2016). Though a dynamite source is relatively simple to model, land surveys typically use vibroseis sources to acquire large volumes of data more quickly, but unfortunately vibroseis source signatures are more difficult to estimate. There have been some recent developments in vibroseis source modelling for multiple sources (Ikelle 2007) and through modelling Green's functions from several locations simultaneously (Neelamani et al. 2008). Although, many methods are approximate and do not fully represent the complex interaction of source with the ground (Rowse & Tinkle 2016), careful calibration of the source has led to successful applications of FWI to vibroseis data (e.g. Plessix et al. 2012). An additional next step that could be useful is to include anisotropy and generate fully anisotropic data, and then investigate whether we can recover the anisotropy (Debens 2015). If we can extract the anisotropic Thomsen parameter ε (Thomsen 1986), we may be able to use this to estimate fracture properties, and even track changes in fracture density and fill over time. Likewise, extending the combined survey to include receivers in the access shaft walls should improve the inversion results and may be valuable should the methods be extended for monitoring (Marelli et al. 2010). In summary, development of methodologies to characterize fracture evolution through combined seismic surveys and 3-D FWI could be powerful in improving our understanding of rock-property evolution relevant to the GDF.

C O N C L U S I O N S
We conclude that 3-D FWI of surface seismic data may be a useful tool in recovering subsurface velocity structure at a potential GDF site, before and after tunnelling. The addition of receivers within the tunnel results in more complete recovery of the EdZ velocity anomaly, and improves the inversion result whether we use reduced offset surveys or include basic or complex tunnels. Notably, 3-D FWI can recover the velocity and shape of the EdZ, features that could not be revealed by other geophysical methods. Site-specific modelling of the GDF and surrounding geology before construction will ensure that the geometry selected for planned seismic acquisition is appropriate, in particular to ensure tunnel receivers are Downloaded from https://academic.oup.com/gji/article-abstract/215/3/2035/5101440 by University of Leeds user on 05 December 2018 placed in any surface survey shadow zones. Importantly, for any future FWI application to seismic data acquired across a GDF, the 3-D codes used here have some additional capabilities that may be important for accurately characterizing the site.

A C K N O W L E D G E M E N T S
We thank Christopher Juhlin and Fengjiao Zhang for their thoughtful reviews which helped improve the manuscript. The research was funded through the Natural Environment Research Council (NERC) Radioactivity and the Environment (RATE) Grant NE/L000423/1 (Hannah Bentham and Doug Angus). We thank GeoRepNet for the Early Career Funding enabling collaborative research visits (Hannah Bentham and Joanna Morgan). We gratefully acknowledge Radioactive Waste Management, Environment Agency and British Geological Survey for guidance and in manuscript preparation. The data and models were processed in ProMAX and Matlab, and figure editing was completed in Inkscape. We thank the sponsors of the FULLWAVE consortium for support in developing the 3-D FWI software used here. Contains data supplied by permission of NERC and University of Leeds hosted by the National Geoscience Data Centre (Bentham et al., 2018).