Interactive comment on “ Electric resistivity and seismic refraction tomography , a challenging joint underwater survey at Äspö Hard Rock Laboratory

We like to thank all the reviewers and the editor for their constructive critics that will help to improve the manuscript significantly. In most cases we have additions to the text, one case (going over several points of reviewer 2) we are going to add another figure about the use of different measures (coverage and resolution radius) for appraising model uncertainty. This will, along with the other changes, help to be a guideline to use geophysics in similar settings and thus be of wide interest for the non-geophysical readership of the journal.

Abstract.Tunnelling below water passages is a challenging task in terms of planning, pre-investigation and construction.Fracture zones in the underlying bedrock lead to low rock quality and thus reduced stability.For natural reasons, they tend to be more frequent at water passages.Ground investigations that provide information on the subsurface are necessary prior to the construction phase, but these can be logistically difficult.Geophysics can help close the gaps between local point information by producing subsurface images.An approach that combines seismic refraction tomography and electrical resistivity tomography has been tested at the Äspö Hard Rock Laboratory (HRL).The aim was to detect fracture zones in a well-known but logistically challenging area from a measuring perspective.
The presented surveys cover a water passage along part of a tunnel that connects surface facilities with an underground test laboratory.The tunnel is approximately 100 m below and 20 m east of the survey line and gives evidence for one major and several minor fracture zones.The geological and general test site conditions, e.g. with strong power line noise from the nearby nuclear power plant, are challenging for geophysical measurements.Co-located positions for seismic and ERT sensors and source positions are used on the 450 m underwater section of the 700 m profile.Because of a large transition zone that appeared in the ERT result and the missing coverage of the seismic data, fracture zones at the southern and northern parts of the underwater passage cannot be detected by separated inversion.Synthetic studies show that significant three-dimensional (3-D) artefacts occur in the ERT model that even exceed the positioning errors of underwater electrodes.The model coverage is closely connected to the resolution and can be used to display the model uncertainty by introducing thresholds to fade-out regions of medium and low resolution.A structural coupling cooperative inversion approach is able to image the northern fracture zone successfully.In addition, previously unknown sedimentary deposits with a significantly large thickness are detected in the otherwise unusually well-documented geological environment.The results significantly improve the imaging of some geologic features, which would have been undetected or misinterpreted otherwise, and combines the images by means of cluster analysis into a conceptual subsurface model.

Introduction
Underground structures have become an increasingly important part of modern infrastructure, and the possibilities to improve construction approaches have attracted much attention.With constantly reduced space for new structures on the surface, underground space is attractive for use in the transportation sector to challenge the increase in traffic in and around cities or as underground storage facilities.Geological uncertainties increase the risk of delays and thus the costs of underground construction.A detailed subsurface model is essential for reducing risks and for a successful project.In order to ensure a smooth construction phase, a critical point is to locate weak zones, especially those that can generate a large inflow of water, causing problems and slowing down the construction progress.Except for southwestern Scania and the islands Gotland and Öland, crystalline bedrock is the dominant material for underground infrastructure construction in Sweden.For these geologic conditions, weakness zones that are im-portant for the underground design are normally indicated by dry, water-bearing or sediment-filled fractures.
Two methods for site investigation in crystalline bedrock are drilling and surface-based or borehole geophysics.Drilling is often the first choice since it provides high resolution and accuracy at any given depth.Nevertheless, drilling is expensive and delivers only point information.Therefore, surface-based geophysical methods have gained more attention, since they provide continuous models that reveal the extreme points and an opportunity for extrapolation into 2-D or 3-D space.The usage of geophysics has increased lately to obtain more continuous and comprehensive subsurface models.Recently, the Swedish transportation authority has provided funding for research in an increasing number of projects with the aim of developing site investigations based on additional geophysical measurements for mapping the structure and quality of the rock mass.Dahlin et al. (1999) report a case in which electrical resistivity tomography (ERT) has been used successfully to map weak and permeable rock in an onshore railway tunnel project in Sweden.Ha et al. (2010) used different geoelectrical applications to detect weak zones of approx.40 m × 40 m during underground construction.In the Norwegian R & D project "Tunnels for the citizens", which was funded by the road administration, several publications (Karlsrud et al., 2003;Palmstrøm et al., 2003;Rønning et al., 2013;Wisén et al., 2012;Lindstrøm and Kveen, 2004) report that elaborate site investigations are important in a controlled tunnelling process, but also that further studies are needed.Rønning et al. (2013) assessed ERT, refraction seismics, very low frequency (VLF) electromagnetics and the AM-AGER method (aeromagnetic and geomorphological relations).They concluded that these are all able to locate fracture zones and state that ERT is able to give more hints to the fracture width, dip and depth extent compared to the other methods used.They also suggest a quantitative rock quality measure on the basis of resistivity values.Refraction seismics has long been an established method for information on fracture width and seismic p-wave velocity; the latter has an obvious coupling to the hardness of the rock and hence to rock quality (Bergman et al., 2006).During an investigation for a road tunnel in Norway, Ganerød et al. (2006) found that 2-D resistivity and refraction seismics are the most suitable geophysical methods, while ERT gave detailed results at lower costs compared to seismics.Diaz et al. (2014) successfully conducted seismic refraction and ERT surveys and associated resistivity and velocity changes with the main and secondary structures of a major fault zone.The final velocity and resistivity models were also consistent with deformed sedimentary units.Another multidisciplinary geophysical approach for mapping a fault zone is given by Malehmir et al. (2016).Heincke et al. (2010) used seismic and electric tomography to assess the rock quality on a hard rock slope in Norway.
Several methods are generally combined to overcome the limits of the natural resolution and corresponding ambiguity in inversion and interpretation.One example for synthetic and field data is given by Garofalo et al. (2015).Seismic data and ERT were used to reduce model ambiguities to improve the estimation of the geophysical parameters.With the joint inversion approach, data sets from different methods are used to constrain each other.The general assumption is that subsurface structures lead to parameter changes for the different methods, i.e. imaging the same underlying geologic conditions.Smoothness constraints often prevent the correct mapping of sharp interfaces.Different joint inversion approaches exist.The cross-gradient method for ERT and seismic data is explained in Gallardo and Meju (2004).Although a significant improvement in the results compared to separated 2-D inversions was observed, it can only be applied on regular grids, which makes an accurate incorporation of topography difficult.Juhojuntii and Kamm (2015) present a joint inversion algorithm using seismic refraction and ERT data, assuming a fixed number of layers.The derived subsurface models agreed well with independent in situ tests, but the authors also stated that their approach would lead to misinterpretations in environments with smooth subsurface variations.
This paper describes a representative case study for the combination of geoelectric and refraction seismics in typical Scandinavian geologic conditions at a coastal region.The survey was conducted at Äspö Hard Rock laboratory (HRL) and designed to perform a joint inversion on the data.The main objective was the localisation and characterisation of fracture zones under challenging conditions, which are the extreme variation in electrode coupling, possible 3-D effects on ERT data and high acoustic damping due to gasbearing sediments.Dahlin and Wisén (2016) and Günther and Südekum (2007) showed that underwater field surveys are possible and quite promising.Loke and Lane (2004) demonstrated that the water layer has a large effect on apparent resistivities, but subsurface resistivity can be recovered if water resistivity and seabed topography are properly incorporated in the finite-element meshes.The methodical approach of a structurally coupled joint inversion presented in this study shows how results can be improved such that an easier and more unique interpretation of the underground models is possible.In order to increase the reliability of the results, a combined inversion and interpretation was investigated.This was done by joint inversion followed by a cluster analysis as an additional integrated interpretation approach.After describing the site conditions and the numerical background, we show a synthetic study on the 3-D effects and the influence of the seabed topography on ERT data before the analysis and interpretation of the field data is presented.

Site description
The Swedish Nuclear Fuel and Waste Management Company (Svensk Kärnbränslehantering AB; SKB) started to design a solution for the deep final disposal of nuclear fuel.Äspö HRL is SKB's underground facility for research that tests concepts for the final disposal of nuclear waste material in hard rock (Rhén et al., 1997).The laboratory has provided a full-scale test environment for different technological solutions.It has now mainly fulfilled its purpose so that the laboratory has also become available for other branches of research.The facility provides a research opportunity in a well-documented and relatively undisturbed geological environment that is representative of many Swedish metropolitan areas.
The Äspö Hard Rock Laboratory is located on the east coast of the Baltic Sea, about 400 km south of Stockholm (see Fig. 1).From 1990 to 1995, the excavation of a 3600 m tunnel that connects the nuclear power plant with the disposal at approximately 450 m of depth was conducted.During the construction phase, a detailed site characterisation was carried out that included geological, hydrogeological and geochemical investigations.
The Äspö bedrock is part of the Transscandinavian Igneous Belt (TIB) that extends from southern Sweden toward the north and northwest.Generally, granitoids and volcanic rocks can be found in the TIB.Four rock types are dominant: the Äspö diorites, Ävrö granite, greenstone and finegrained granite.Wikberg et al. (1991) found that continuous magma mixing processes supported the development of dikes and mafic inclusions, which form an inhomogeneous rock mass.The crystalline bedrock exhibits porosities of 0.4-0.45% for the Äspö diorite and 0.23-0.27% for the fine-grained granite (Stanfors et al., 1999).During the preinvestigation of Äspö HRL, fracture zones were divided into major (width > 5 m) and minor (width < 5 m) categories.The majority of the fractures are oriented northwest-southeast (Berglund et al., 2003).All fracture zones that are important for this field survey are depicted as black lines in Fig. 1.
The filling material of the fractures was extracted from drill cores and analysed.Missing unconsolidated material that might have been additionally filling the fractures was probably washed away and thus not taken into account in these analyses.The crystallised calcite in the fractures was possibly formed by hydrothermal processes and can be used as an indicator for water paths in the rock (Wikberg et al., 1991).This indicated that fractures in the N-S and E-W directions most likely conduct or formerly conducted water.According to Wikberg et al. (1991) all fracture zones are at least partly water bearing.They also gave a judgement of the fracture zones according to Bäckblom et al. (1990).Based on that, the most critical fracture zone along the measured profile is NE-1, which is judged as "certain".EW-3 is also judged "certain", but hydraulically of minor importance.NE-3 and NE-4 are judged as "certain" as well.Both consist of several subzones that are one to a few metres wide, some of which are open fractures that are hydraulically highly conductive.In general, the fracture zones NE-3, NE-4 and EW-7 are judged to be "probable" in a hydraulic sense (Wikberg et al., 1991).The authors also stated that the Quaternary sediments on top of the bedrock were supposed to be scarce at the Äspö test site.Due to the deep target of the Äspö HRL within the bedrock, no detailed investigation of the Quaternary sediments was carried out.Vidstrand (2003) stated that the unconsolidated overburden should rarely exceed 5 m in thickness and consists mainly of clay, sand and gravel.

Electrical resistivity tomography
ERT measurements were carried out along a profile in the N-S direction simultaneously with the seismic survey on 20-24 April 2015.The profile lies between Hålö and Äspö (see Fig. 1) to the west of the tunnel line, about 10 m away from a small island.Electrodes were placed onshore and underwater with a 5 m electrode spacing along a 780 m profile.Data were recorded using the multichannel instrument ABEM Terrameter LS (Guideline Geo, Sundbyberg, Sweden).A multiple-gradient array (Dahlin and Zhou, 2006) was employed to ensure fast measuring progress as it can fully exploit the recording channels.The resistivity of the water was measured with a micro Wenner alpha array at different depths with the ABEM Terrameter LS.The collected ERT data were first published in Fennvik (2015).A model based on accurate bathymetry measurements was used to determine the heights of the sensor positions at the seabed.A nearby power plant caused a high noise level in the ERT data.Large variations in the contact impedance between the water www.solid-earth.net/8/671/2017/Solid Earth, 8, 671-682, 2017 and the rock outcrops created a technically difficult measuring situation.Contact resistances, including cable resistance, started from 100 for electrodes in brackish water and exceeded 100 k on rock outcrops.An over-amplification of the signals was avoided due to an automated gain control of the instrument.Furthermore, the input channels are galvanically separated; i.e. one channel can have a high gain and the next channel a low gain, avoiding any problems.
The full wave form of the transmitted and received signals was recorded in order to recover possibly valuable IP signals from the data.However, the signal-to-noise ratio was sufficiently good for recovering DC resistivity but not IP data.About 6700 data points were gathered during the ERT survey.While the raw data were processed, combinations with uncoupled electrodes were identified and all combinations containing these electrodes were deleted.To account for the variable data quality of the individual data, usually a data error is estimated by a fixed percentage and a voltage error.They can be retrieved by analysing reciprocal measurements (Udphuay et al., 2011), which were not available here.Therefore, we used the default values of 3 % noise and a voltage error of 0.1 mV.

Seismic refraction tomography
The green dashed line in Fig. 1 marks the profile for the seismic refraction.Hydrophone streamers were laid out with 91 hydrophones in total and a 5 m spacing along a 450 m profile line.For data acquisition, the instruments ABEM Terraloc (Guideline Geo) and Geometrics Stratavisor (San Jose, CA, USA) were used, both with 48 channels and a 5-channel overlap of the two streamers.Hydrophone positions were determined by a differential GNSS, while the topography of the seabed was mapped with a multibeam echo sounder (MBES).For all underwater sensors (electrodes and hydrophones), a very accurate DTM (digital terrain model) from the MBES survey was used for the heights.The positions of the sensors (E/S) are coincident and measured with sufficient accuracy.For the excitation of seismic p-waves, small explosives were placed approximately 0.5 m above the seabed.Shots were performed every 20 m.Data were first processed and published in Lasheras Maas (2015).Due to time constraints, not all planned shots were fired, and hence there are two small gaps in the data coverage in the northern part of the data set.Raw data processing revealed that the seismic signal quality was significantly reduced in the southern part of the profile, which made it difficult to pick first arrivals.However, no additional filters were used during the raw data processing.About 650 first-arrival times were semi-automatically picked and manually checked using the software package Rayfract (http://www.rayfract.com).

Numerical modelling and inversion
We used the open-source ERT software package BERT (Boundless Electrical Resistivity Tomography) for ERT inversion (Günther et al., 2006b) using irregular triangle meshes to accurately take into account both the surface and submarine topography (Rücker et al., 2006).Furthermore, we used the underlying framework pyGIMLi (Python Geophysical Inversion and Modelling Library; http://www.pygimli.org)for the refraction tomography and the implementation of the coupled inversion.

Inversion
Geophysical inversion describes the process of estimating a model with a forward response that fits the observed data.
The linearised problem for ERT is given in Eq. ( 1a) and for seismic in Eq. ( 1b).Here, the model parameters are either the logarithmic resistivities or the velocity/slowness held in the model vector m: The Jacobian matrices J and A contain the partial derivatives ∂ρ a,i /∂m j (ERT) or ∂t i /∂m j .Apparent resistivities (ρ a ) are held in d and travel times in t.The inversion of ERT and SRT (seismic refraction tomography) was done by a smoothnessconstrained minimisation using the cost function containing an error-weighted data misfit d and a model roughness m weighted by the regularisation parameter λ.
As the travel time t between source and receiver along a ray path is given by t = n i=1 l i /v i , it is a linear combination of the path length l i and the slowness 1/v i for a segment i.The difference between the individual data points d i and the corresponding forward responses f i (m), both as logarithmic apparent resistivities or travel times, is weighted by their individual errors i .
The roughness (second term in Eq. 2b) consists of the derivative matrix C applied to the model m (Günther et al., 2006b), whereas each row in C is associated with a boundary.Additional model constraints can be incorporated in the object function by extending m from Eq. ( 2b) with a weighted model functional W c (Rücker, 2011), resulting in The weighting matrix W c is diagonal and contains the elements w i , representing penalty factors for the different model cell boundaries (Günther et al., 2006a).Very small values can Cluster analysis

Resistivity Refraction
Figure 2. Scheme of the coupled inversion approach, where the roughness C of one inversion is influenced by the other (Günther et al., 2006a).
lead to sharper boundaries.The limited amount and quality of recorded data leads to a non-unique inversion result.Due to the model smoothing needed for mixed determined problems, it is possible that sharp boundaries appear as transition zones that lead to misinterpretations.A structurally coupled joint inversion finds common structures and allows the models to emphasise these and reduce smoothing effects (Gallardo and Meju, 2004).Here, the roughness vector r = CW m m is used to calculate the mutual penalty factors w i using Eq. ( 4) after Günther et al. (2010): The parameters a, b and c are used to adjust the coupling strength and the influence of the gradients.Differently from the latter approach, we multiply the w i of the different methods and calculate one weighting matrix for both methods.
A certain number of separated iterations is done before the coupling starts so that each method can first independently develop structures before their similarity is promoted.A schematic sketch of the structurally coupled joint inversion is shown in Fig. 2.
Forward modelling and inversion are done accurately on unstructured finite-element (FE) meshes that allow for the incorporation of both the surface and the underwater topography.The finite-element mesh used for the joint inversion of seismic and ERT data is shown in Fig. 3.
The shown mesh consists of three regions that present the background (red), the water (blue) and the parameter domain (green) on which the data inversion is conducted.The orig- inal mesh extension is 1250 m in the x and approximately 420 m in the z direction and is clipped for display reasons.
In situ water conductivity measurements showed resistivity values of about 1.4 m and negligible variation with position or depth.As the seismic velocity of water is constant (about 1400 m s −1 ), the water region can be assumed as homogeneous and is incorporated as a single region with a fixed resistivity or velocity so that the correct values are used for the forward calculation but are not subject to inversion.The parameter domain is extended to approximately 790 m in the x and 190 m in the z direction.Additionally, an outer background region is needed for accurate forward calculation using approximate boundary conditions (Rücker et al., 2006).Although the seismic line is shorter than the ERT, the shown mesh was used for both data sets in the joint inversion.The parameter domain consists of about 3500 cells, which is the number of model parameters.More details on region-based inversion can be found in Rücker (2011).

Model appraisal and display
All shown inversion results are faded out using the coverage to point out the contribution of the model parts to the data.The calculation of the coverage is based on the sensitivity, which is the partial derivative S i,j (m n ) = ∂f i (m n ) ∂m j .Whereas m = log ρ are the model parameters and f = log ρ a is the forward response (both logarithmically transformed) for ERT, the seismic model parameters are m = 1/v (slowness) with the corresponding forward response f = t (travel times).The summation of all sensitivities for each model parameter gives the coverage for the model cell assigned with this parameter.Unlike ERT, a normalised coverage is calculated for seismics, which is either 1 or 0 depending on whether a ray crosses a cell or not.Additionally, resolution radii after Friedel (2003) were calculated for the purpose of a comparison with the coverage.For that the model resolution matrix R M is required, which can be calculated by a singular value decomposition (SVD) of the Jacobian matrix of the final model, i.e. the resistivity or velocity/slowness distribution.The radius r j = A cell /(R M j,j π ) is calculated for each model cell j and is an equivalent of the cell area A cell to a  In general, low-resolution radii correspond with a high coverage and vice versa.The water and the low-resistive sediments between x = 180 m and x = 550 m lead to a reduced investigation depth with high-resolution radii starting at approx.50 m of depth.Compared to that, the onshore model parts at x < 180 m and x > 550 m show a medium reliability up to 70 m of depth with resolution radii < 50 m.Figure 4 clearly shows the existing relationship between the coverage and the model resolution radii.Thus, the coverage can also be used as a resolution measure to display well-, mediumand poorly resolved model regions.In the following, two coverage thresholds were used to define regions of high, medium and low certainty.The low-certain region is completely blanked out and considered as untrustworthy, while the region of high certainty is imaged without shading and can thus be judged as trustworthy.

Synthetic study on 3-D effects and seabed topography
We follow a strict two-dimensional (2-D) scheme, i.e. assuming constant values perpendicular to the profile.For the given test site, it can be assumed that the seismic refraction data are not or only minimal corrupted by 3-D effects.By picking first arrivals, only signals that took the shortest way or travelled in the fastest medium are taken into account.The small island next to the profile consists of the same bedrock as that directly in the profile line.Assuming the same velocity, the recorded first arrival is still from the signal travelling in the profile line, because it is the shortest way.The small bay (water body) north of the profile would be a low-velocity anomaly because the velocity in water is lower than in bedrock.Therefore, the bay can be ignored because a refraction only appears for an increasing velocity.Three-dimensional (3-D) effects occur if significant resistivity changes perpendicular to a 2-D profile are present.According to the test site map in Fig. 1, severe 3-D effects can be expected near the small island in the middle of the profile and in the northern part, where the water continues just a few metres next to the profile.The latter is not expected to have a significant effect on the first-arrival times, since these are related to the smallest distance to the layers.It will, however, have an effect on the measured apparent resistivity by all materials present within the measured volume.In order to appraise the expected shapes and magnitudes of 3-D effects, we generated a simplified model based on the Äspö geometry.The underlying model used for generating synthetic data is shown in Fig. 5.The water body is simulated by a cube with an extension of 450 m in the x direction starting at x = 100 m, being 10 m in depth and infinite in the y direction.
A large cube simulating the bedrock (brown) surrounds the water cube, with an infinite extension in the x, y and z direction.The water (blue) is assigned with a resistivity of 3 m, while the bedrock is assigned with 3000 m.Two anomalies are inserted representing the island in the middle and the small bay at the northern end of the ERT profile.The island (red) is a 10 m thick cube, with an extension of 90 m in the x and 70 m in the y direction, placed between x = 370 m and x = 460 m with a distance of 10 m to the ERT profile.The small bay at the northern part (green) is incorporated by a rectangular cube with an edge length of 100 m (x, y direction) and 3 m of depth.It starts directly after the water cube at x = 550 m with a distance of 5 m to the profile.The ERT line consists of 153 electrodes, starts at x = 15 m and y = 0 m and is aligned along the x direction.The simulated survey is identical to the field measurements except that the electrodes are assumed to be at the surface and topography is neglected.
The ERT profile is marked with red spheres in Fig. 5.For reference, we additionally calculated data from a 2-D model, where the island is assigned with a water resistivity of 3 m and the bay with 3000 m (bedrock), i.e. with no 3-D effects.Both data sets were corrupted with Gaussian noise with an error level consisting of 3 % plus a voltage error of 100 µV.A smoothness-constrained inversion was performed to estimate resistivity models from the two synthetic data sets.Figure 6a and b show the inversion results from the data set with and without 3-D effects.The ratio between those two is shown in Fig. 6c.
Figure 6a shows the expected smooth resistivity distribution with a horizontal interface between the simulated bedrock and water.When the island and the small bay are included in the underlying model, serious 3-D effects occur.Äspö case, a water-filled valley with a depth of 10 m and a length of 550 m was used.The reference model contains a flat seabed, whereas cases one and two contain a depth variation of ±0.30 m.For the first case, the depth of the seabed was set to −10.3 m between x = 230 m and x = 300 m and to −9.7 m between x = 380 m and x = 450 m.The depth varies, alternating from −10.3 to 9.7 m for 250 m ≤ x ≤ 395 m for the second case.While the data were generated, 2 % Gaussian noise was added.Afterwards, geometric factors for the first and second case were replaced with the reference data in order to simulate a flat seabed.One mesh for the data inversion was used with a flat seabed.The ratios between the two cases and the reference data set are shown as pseudo-sections of the simulated data in Fig. 7 The deviation due to the changed seabed topography is in the range of ±20 %. Figure 7a shows a clear pattern due to the changed model, which is a lower apparent resistivity for a slight downwards shift of the seabed and an increased apparent resistivity for an upwards shift.Compared to that, an alternately varying seabed leads to a rather random pattern (Fig. 7b).This simple synthetic study confirms that the ERT data gathered at the Äspö test site are contaminated or distorted by 3-D effects that have to be taken into account when interpreting the results.
www.solid-earth.net/8/671/2017/Solid Earth, 8, 671-682, 2017 A smoothness-constrained inversion was done with the abort criterion χ 2 = d /N = 1; i.e. the data are fitted within their errors.A visual inspection of the data misfit ensured that there was no more unresolved structure.The L 1 norm data (robust) inversion was used to account for remaining outliers in the ERT data set that lead to poor data fits.Nevertheless, the apparent resistivities cover several orders of magnitude (3-47 000 m) and extraordinarily high resistivity variations occur, which is challenging for ERT inversion.The ERT inversion result is shown in Fig. 8a using the coverage (sum of the absolute Jacobian values over all data for each model) for alpha shading.In the middle of the profile, the penetration depth is limited due to the low-resistive water body and the anomalies below.
Outcrops of the bedrock lead to high resistivities of about 35 000 m at the northern and southern ends of the profile.A low-resistive zone appears at x = 200-600 m directly below the sea.The depth varies between approximately 80 m at x = 270 m and 30 m at x = 450-600 m.As such a deep weathering zone seems implausible and the resistivity is too low for usual weathering, we interpret this structure as a deep valley filled with sediments.This has not been documented by previous investigations conducted in the construction phase of the test nuclear waste disposal.The lowresistive zone is extended diagonally downwards towards the north for x > 600 m at a depth range of 50-100 m.Although the coverage is low for this part, it is still possible that this feature indicates fractured water-bearing bedrock.
Resistivities of about 500 m at x = 100-200 m and a depth of 100 m indicate a larger transition zone that continues below the sediment body.This could possibly lead to an incorrect depth of the sediment-filled valley and thus the bedrock interface.It also prevents any further interpretation regarding possible fracture zones.
The inversion result of the refraction seismics shown in Fig. 8b images the interface to the bedrock more accurately.However, the poor signal quality in the southern part results in a lower coverage and thus larger uncertainty.To display the inversion result, a standardised coverage was calculated, which is either 0 or 1 depending on whether any ray travels through a model cell or not.
According to Fig. 8, the crystalline bedrock appears as a high-velocity zone of about 5600 m s −1 , which agrees with the velocity for intact crystalline rock at Äspö HRL given by Wikberg et al. (1991). Brodic et al. (2016) recently showed that the velocity decreases from > 5000 m s −1 for intact rock down to approx.4200-4700 m s −1 for fracture zones.sediments exhibit a minimum velocity of about 1000 m s −1 , which is below the velocity of water (1400 m s −1 ).A reason could be gas contained in the sediments, which reduces the acoustic velocity for frequencies below 1 kHz (Wilkens and Richardson, 1998).This is supported by the presence of gas bubbles rising to the water surface during the blasting.Gasbearing sediments were also reported by Dahlin et al. (2014) near Stockholm, which has a similar geologic history.It is assumed that the gas-bearing sediments lead to the poor data quality in the southern part by damping the seismic signals.
No further low-velocity zones appear at larger depths.
To summarise, a (possibly gas-bearing) sediment body could be identified, which appears as a zone of low resistivities and velocities.Furthermore, the interface towards the bedrock could be found by the joint interpretation of the separated inversion results.However, the bedrock appears with a low resistivity due to the large transition zone.Fracture zones are not visible in the separated inversion results (Fig. 8) because of a low coverage in the refraction model and a large transition zone in the resistivity model.
In order to improve the results and enable further interpretation, a structurally coupled joint inversion of the ERT and seismic data was performed.To ensure that common structures are present in the models, the first four iterations were done separately.A robust data fit, i.e.L 1 norm, was used for ERT data inversion, while the first arrivals were fitted using the L 2 norm (least squares).Both data sets were fitted within their errors, i.e. with χ 2 = 1.1 for ERT and χ 2 = 1.3 for refraction data.In this case, the RMSE (root mean square error) for the first-arrival fit was about 2.4 ms.The result is shown in Fig. 9.Both models show significant changes compared to the separated inversions and allow for further interpretations.Generally, most changes occur in the resistivity model, while the velocity model shows only small improvements.The lowresistive zone, which corresponds to the sediment-filled valley, appears thinner followed by a much smaller transition zone.This reduces the ambiguity in estimating the bedrock interface.The bedrock is also assigned with a higher resistivity, which is more realistic as it agrees with the resistivity of near-surface rock outcrops at the northern and southern ends of the profile.
Additional structural constraints that moved from the velocity to the resistivity model pointed out the diagonal lowresistive zone in the northern part in more detail.This anomaly matches very well with the water-bearing fracture zone NE-1 in the northern part of the profile.The southern fracture zones NE-3, NE-4 and EW-1 cannot be identified directly.Possible explanations could be that these are (i) too small to be detected from the surface or (ii) filled with a material so that no parameter contrast appears.
According to the synthetic study, the low-resistive feature directly at the surface at x = 610 m and, in part, the diagonal low-resistive zone at x = 600 m are most likely caused by 3-D effects and should not be interpreted any further.Following Günther et al. (2006a), a post-processing of the two inversion models was done using a cluster analysis to obtain a simplified result (Fig. 10).For clustering the resistivity and velocity model, a modified mean-shift algorithm approach was used, which is described in Comaniciu and Meer (2002).The input for this algorithm is a feature space that consists in this case of resistivities and velocities.In order to analyse the feature space, a window or bandwidth is needed.The bandwidth can be determined by an estimator that uses a selected quantile as input, whereas the quantile is defined between 0 and 1.In general, a low quantile will produce a larger number of clusters than a high quantile.In contrast to cluster number-driven algorithms such as the K-means algorithm (see Joydeep and Alexander, 2009), the input is data and a window to the data.Therefore, the selection of clusters is driven only by data and not by an arbitrary number of clusters.
As data input for the clustering, we only used model parameters included by the coverage of the seismic result (displayed cells in Fig. 9b) because the seismically covered volume is also covered by ERT.
The data-driven cluster algorithm divided the model parameters into three clusters that represent sedimentary deposits, the bedrock and the transition zone between those two.It can most likely be assumed that the interface between the sediments and the bedrock is within the third cluster.
www.solid-earth.net/8/671/2017/Solid Earth, 8, 671-682, 2017 As a final interpretation of the presented ERT and seismic results, a conceptual model was developed (Fig. 11).The primary origin of the deep sedimentary deposits can be explained by glacial erosion.The small valley was formed between the fracture zones NE-3 and NE-4.It might have been easier to erode the bedrock along zones with an already low rock quality.Two possible explanations can be given for the remaining transition zone at the bottom of the sedimentary valley.The first is that the bedrock-sediment interface is (i) fractured or weathered to a certain extent, and the second is (ii) that coarse sediments could have been deposited before fine-grained marine material was sedimented above.The latter possibility is visualised by the dark yellow and orange parts at the bottom of the valley in Fig. 11.As the medium velocities north of the sedimentary valley appear slightly thicker, the most probable explanation could be weathered bedrock.During an earlier investigation, it was found that the NE-1 fracture zone in the northern part of the model is water bearing at its boundaries and dry in its core due to clay deposits.Thus, it appears as a zone of lower resistivities and velocities.Only the NE-1 fracture zone could be identified by this survey, although the fracture zones NE-3, NE-4 and EW-3 are also partly water bearing according to Wikberg et al. (1991).As shown in Fig. 1, NE-4 and EW-7 are close to each other at the profile line, which means that they most likely cannot be imaged separately by ERT measurements.In addition, the low-resistive sediments are a complicating factor that may mask the fracture zones by reducing the model resolution such that it is not sufficient to resolve the fracture zones NE-3, NE-4 and EW-7.

Conclusions and outlook
A combination of refraction seismics and ERT data has been tested on an underwater profile crossing a water passage along part of the access tunnel that connects surface facilities with an underground test laboratory at the Äspö Hard Rock Laboratory.The aim was to detect fracture zones in a well-known but logistically challenging area.Co-located sensor positions for ERT and seismics were used on a 450 m underwater section of the 700 m ERT profile.
A synthetic study inspired by the geologic conditions of the Äspö test site showed that significant 3-D effects are expected that contaminate the ERT data and thus influence the obtained inversion result.This was taken into account to prevent the misinterpretation of the final inversion results.The results of the separated inversions showed a previously unknown sediment-filled valley that appeared as a zone with low resistivities and low velocities, even in an unusually well-documented geological environment.The poor coverage of the seismic model in the northern and southern parts of the profile in conjunction with the large transition zone of the ERT result prevent further detailed interpretations.However, the water-bearing fracture zone NE-1 could be iden-tified by the results of the structurally coupled joint inversion.The evaluation shows that the joint inversion approach combining ERT and seismics has very promising results for three reasons: (i) the decreased extent of the transition zone, (ii) the more reliable interpretation of two independent parameters and (iii) their combination by a clustering approach.Although the refraction seismic does not cover the fraction zone NE-1, the additional constraints by the joint inversion helped to determine the fracture zone with ERT.The southern fracture zones NE-3, NE-4 and EW-1 could not be detected due to the missing parameter contrast and/or the model resolution.The latter is considered to be the main reason, which was shown by the distribution of the model resolution radii and coverage.The reduced investigation depth of ERT is due to the fact that the current preferably flows through low-resistive bodies (water or sediments) and is the major disadvantage of this method.
The comparison of the joint inversion with the separated inversion result shows significant improvements.Therefore, the combination of geoelectric and seismic refraction is recommended as the standard tool for site investigations under geologic conditions similar to those presented.

Figure 1 .
Figure 1.The location, major fracture zones (black lines) after Stanfors et al. (1999), the scheduled ERT profile (solid red line) and the seismic (dashed green) line at Äspö Hard Rock Laboratory.

Figure 3 .
Figure3.Cropped mesh used for the structurally coupled inversion of ERT and seismic data.Three regions are used: (i) the background (in red; much bigger) to prevent influences of the boundaries, (ii) the parameter domain (green) on which the inversion is done and (iii) the water region (blue), which was fixed.

Figure 4 .
Figure 4. Calculated model resolution radii (a) and coverage (b) for ERT using the Jacobian matrix of the final model.

Figure 5 .Figure 6 .Figure 7 .
Figure 5. Sketch of the synthetic model used to generate synthetic data.It reflects a simplified version of the Äspö test site conditions.The red spheres mark electrode positions, the blue coloured areas simulate a low-resistive body, like sea water, and the brown parts mark highly resistive bodies, like bedrock.

Figure 8 .
Figure 8. Separated inversion results of the ERT data set (a) and the refraction seismic data (b).The shading is based on the coverage.

Figure 9 .
Figure 9. Joint inversion result with resistivity (top panel) and velocity (bottom panel) distribution.The shading is based on the coverage of each model cell.

]Figure 10 .Figure 11 .
Figure 10.Cluster analysis of the joint inversion result using tree clusters.The upper picture shows the spatial distribution of the clusters and the lower one shows the parameter distribution within each cluster.