# Joint inversion of lake-floor electrical resistivity tomography and boat-towed radio-magnetotelluric data illustrated on synthetic data and an application from the Äspö Hard Rock Laboratory site, Sweden

Joint inversion of lake-floor electrical resistivity tomography and boat-towed... Summary The electrical resistivity tomography (ERT) method provides moderately good constraints for both conductive and resistive structures, while the radio-magnetotelluric (RMT) method is well suited to constrain conductive structures. Additionally, RMT and ERT data may have different target coverage and are differently affected by various types of noise. Hence, joint inversion of RMT and ERT data sets may provide a better constrained model as compared to individual inversions. In this study, joint inversion of boat-towed RMT and lake-floor ERT data has for the first time been formulated and implemented. The implementation was tested on both synthetic and field data sets incorporating RMT transverse electrical mode and ERT data. Results from synthetic data demonstrate that the joint inversion yields models with better resolution compared with individual inversions. A case study from an area adjacent to the Äspö Hard Rock Laboratory (HRL) in southeastern Sweden was used to demonstrate the implementation of the method. A 790-m-long profile comprising lake-floor ERT and boat-towed RMT data combined with partial land data was used for this purpose. Joint inversions with and without weighting (applied to different data sets, vertical and horizontal model smoothness) as well as constrained joint inversions incorporating bathymetry data and water resistivity measurements were performed. The resulting models delineate subsurface structures such as a major northeasterly directed fracture system, which is observed in the HRL facility underground and confirmed by boreholes. A previously uncertain weakness zone, likely a fracture system in the northern part of the profile, is inferred in this study. The fractures are highly saturated with saline water, which make them good targets of resistivity-based geophysical methods. Nevertheless, conductive sediments overlain by the lake water add further difficulty to resolve these deep fracture zones. Therefore, the joint inversion of RMT and ERT data particularly helps to improve the resolution of the resistivity models in areas where the profile traverses shallow water and land sections. Our modification of the joint inversion of RMT and ERT data improves the study of geological units underneath shallow water bodies where underground infrastructures are planned. Thus, it allows better planning and mitigating the risks and costs associated with conductive weakness zones. Electrical resistivity tomography (ERT), Radio-magnetotellurics (RMT), Joint inversion, Fractures, faults, and high strain deformation zones 1 INTRODUCTION Due to the non-linearity and ill-posedness of geophysical inverse problems, seeking a unique model based on geophysical data is a difficult task. Joint inversion is a useful technique to narrow down the range of possible models and has become increasingly popular for interpreting geophysical data sets related to both common and disparate material properties. Two strategies, direct parameter coupling and structural coupling, are mainly used for joint inversion with different types of geophysical model parameters (e.g. Moorkamp et al.2011; Abubakar et al.2012; Haber & Gazit 2013; Moorkamp et al.2016). Direct parameter coupled joint inversion is based on petrophysical parameter measurements in the laboratory and only valid in limited parts of the model domain (Lines et al.1988; Coutant et al.2012). Structurally coupled joint inversion using cross-gradient constraints (Gallardo & Meju, 2003, 2004, 2007, 2011; Tryggvason & Linde 2006; Hu et al.2009; Moorkamp et al.2010; Moorkamp et al.2016) is easy to implement and offers sufficient flexibility to reduce structural coupling where this is in contradiction to the data. Single-property joint inversion is naturally coupled by inverting data from different methods for the same type of material property, for example among different types of electromagnetic (EM) methods (Meju 1996; Kalscheuer et al.2015). For the joint inversion of radio-magnetotelluric (RMT) and electrical resistivity tomography (ERT) data, many relevant publications exist on algorithm developments and case studies. First, Vozoff & Jupp (1975) described the advantages of joint inversion of MT and ERT data by synthetic and field examples using 1-D models. Sasaki (1989) implemented a 2-D joint inversion of MT and ERT dipole–dipole data and also tested it on synthetic and field data sets. Monteiro Santos et al. (2006) evaluated a geothermal field using 2-D joint inversion of MT and ERT data by the technique proposed by Sasaki (1989). Candansayar & Tezkan (2008) developed an algorithm for 2-D joint inversion of RMT and ERT data, which improves the resolution in modeling a fracture zone both in synthetic and field cases. This algorithm was also used to evaluate the impact of sewage irrigation and groundwater contamination in India by Yogeshwar et al. (2012). Kalscheuer et al. (2010) analysed the resolution and variance properties of 2-D resistivity models derived from single and joint inversions of ERT and RMT measurements. Bastani et al. (2012) investigated an oil contamination field in Italy by joint 2-D inversion of ERT and RMT data and compared the results with 3-D single inversion of RMT data. By studying a quick-clay site in Norway, Kalscheuer et al. (2013) demonstrated that the detectability of quick-clay zones could be improved by jointly inverting ERT and RMT data. For a few reasons, joint inversion of RMT (also MT data in general) and ERT data sets can lead to inverse models that are better constrained than those from the inversion of individual data sets. Due to spatial limitations, the ERT measurements may only constrain the shallow part of the subsurface (upper 20–30 m depending on the case) using electrode arrays that are a couple of hundred metres long. In cases with longer electrode arrays (>1 km) and large electrode spacing, ERT penetration depth may be substantially larger, but resolution for the shallow subsurface may be reduced. Hence, RMT measurements with a fixed range of signal frequencies of 14–250 kHz may constrain either correspondingly deeper or shallower parts of the subsurface. Thus, combining ERT and RMT data sets can lead to improved depth coverage of geological targets (Candansayar & Tezkan 2008). Furthermore, the sensitivities of the different methods complement each other; for example, ERT can resolve both resistive and conductive structures moderately well, while RMT is superior at resolving conductive structures but has negligible sensitivity to the resistivity of thin resistive units (Vozoff & Jupp 1975; Candansayar & Tezkan 2008; Kalscheuer et al.2010). Also, measurements by the different methods are differently affected by various types of field noise and couple differently to subsurface structures in the survey area. Finally, joint inversion of RMT and ERT data may to some extent help to determine electrical anisotropy parameters (Christensen 1998). Owing to the improvements in model constraints offered by joint inversions of ERT and RMT data and the increasing number of geo-engineering applications to investigate rock mass underneath shallow water bodies, we extended the joint inversion algorithm by Kalscheuer et al. (2010) to invert boat-towed RMT and lake-floor ERT data. The results of synthetic tests and a real-field application are presented in this work. This is a unique implementation and application, although lake-floor ERT (Dahlin et al.2014; Loke & Lane 2004) and boat-towed RMT (Bastani et al.2015; Mehta et al.2017) experiments have separately been shown, and joint inversion of ERT and RMT data collected on land has been done before (Candansayar & Tezkan 2008; Yogeshwar et al.2012; Bastani et al.2012; Kalscheuer et al.2013). To further improve the estimated joint inversion model, weighting techniques for the different data sets and incorporation of bathymetric data and water resistivity measurements as a priori constraints were employed. We applied the method to a real data set collected in an experiment at the Äspö Hard Rock Laboratory (HRL) 30 km north of the city of Oskarshamn in southeastern Sweden. The target was mainly a dipping fracture system (consisting of several subsets) that was intersected by the access tunnel of the HRL facility and was speculated to reach the lake floor: no clear evidence for this hypothesis was available from previous geoscientific studies. 2 THEORY 2.1 Radio-magnetotelluric method RMT is one type of passive-source EM methods where the signal sources are distant radio transmitters operating in the frequency range of 14 kHz –1 MHz. At sufficiently long distances, the EM signals are considered as plane waves simplifying the estimation of the electrical resistivity of near-surface structures (Bastani 2001; Bastani et al.2012; Tezkan et al.1996, 2000; Turberg et al.1994). RMT data comprise measurements of three components of the magnetic field (Hx, Hy, Hz) and two horizontal components of the electric field (Ex, Ey). In the frequency domain, the electric and magnetic field components are related through the complex-valued and frequency-dependent impedance tensor Z given as   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{E_x}}\\ {{E_y}} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{Z_{xx}}}&{{Z_{xy}}}\\ {{Z_{yx}}}&{{Z_{yy}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{H_x}}\\ {{H_y}} \end{array}} \right],$$ (1)and the vertical magnetic field and horizontal magnetic fields are related through the complex-valued and frequency-dependent vertical magnetic transfer function T given as   $${H_z} = \left[ {\begin{array}{*{20}{c}} {{T_x}}&{{T_y}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{H_x}}\\ {{H_y}} \end{array}} \right],$$ (2)where x and y represent measurement directions in the Cartesian coordinate system. The skin depth   $${\delta = 503\sqrt {\rho /f} },$$ (3)where f is frequency (unit is Hz) and ρ is resistivity (unit is Ωm), is typically used to evaluate the maximal exploration depth at the given signal frequency as 1.5 δ (Spies, 1989). In a layered medium, the resistivity ρ is replaced by an effective resistivity $$\tilde{\rho }$$ (Spies 1989; Huang 2005). In 2-D media, the diagonal impedance tensor elements are zero (Zxx = Zyy = 0) given that x is the strike direction of the 2-D geological structure and y is the profile direction. The apparent resistivities (ρxy/yx) and phases (φxy/yx) defined by Zxy (transverse electrical, TE mode, with currents flowing in x-direction) and Zyx (transverse magnetic, TM mode, with currents flowing in y- and z-directions) are used to estimate the general resistivity distribution:   $${{\rho _{xy/yx}} = \frac{1}{{\mu_0 \omega }}}{\left| {{Z_{xy/yx}}} \right|^2},$$ (4)  $${\varphi _{xy/yx}} = {\tan ^{ - 1}}\left( {\frac{{{\mathop{\rm Im}\nolimits} \left( {{Z_{xy/yx}}} \right)}}{{{\mathop{\rm Re}\nolimits} \left( {{Z_{xy/yx}}} \right)}}} \right)$$ (5)where μ0 and ω are permeability of free air and angular frequency. In areas where the strike direction is not well defined, inversion of determinant impedance data is more suitable to suppress 3-D effects in 2-D models (Perdersen & Engels 2005) than inversion of biased TE- or TM-mode data. Here, the determinant impedance is defined as   $${Z_{{\rm{DET}}}} = \sqrt {{Z_{xx}}{Z_{yy}} - {Z_{xy}}{Z_{yx}}} .$$ (6) 2.2 Electrical resistivity tomography Use of ERT dates back as far as the beginning of geophysics and originally was used to distinguish oil-saturated from water-saturated formations in hydrocarbon exploration (Schlumberger 1939). In near-surface geophysics, the ERT method has been used for hydrogeological, mining and geohazard investigations (Daily et al.2005). In ERT measurements, a direct current (I) is injected into the ground by means of two electrodes (one can be far away from the measurement area), and a potential difference or voltage (U) owing to this current flow is measured between two potential electrodes. The Wenner, dipole–dipole, Schlumberger and gradient (Dahlin & Zhou 2006) electrode configurations are commonly used for ERT surveys (Zonge et al.2005). ERT data are presented as pseudo-sections of apparent resistivity or resistance, which are defined as   $${\rho _a} = k\frac{U}{I}$$ (7)or   $$R = \frac{U}{I},$$ (8)respectively. Here, k is a geometric factor, which depends on the array type and electrode spacing. These measurements are geometrically weighted averages of the true resistivity distribution of the subsurface. For current injection in a homogeneous Earth, the section of the earth above the median depth of investigation (Loke 1999) has the same influence on the measured potential as the section below. This determines roughly to which depth anomalous geological structures can be investigated using a given array type (Loke 1999), and the maximal depth of investigation differs from array to array. 2.3 Inverse theory We used the code Electro-Magnetic Inversion using Least Intricate Algorithms (EMILIA, Kalscheuer et al.2008, 2010) for our implementation, carrying out synthetic tests and inverting the field data. EMILIA uses the finite-difference method to solve the Poisson equation (Dey & Morrison 1979; Kalscheuer et al.2010) and the Helmholtz equations (Siripunvaraporn & Egbert 2000; Kalscheuer et al.2008) for the ERT and RMT methods, respectively. In an inverse problem, an objective function of the form (Menke 1989; Kalscheuer et al.2010):   \begin{eqnarray} \Phi \left( {{{\boldsymbol m}},{\lambda }} \right) &=& {\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{\boldsymbol m}} \right]} \right)^{T}}{{\boldsymbol W}}_{d}^{T}{{{\boldsymbol W}}_{d}}\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{\boldsymbol m}} \right]} \right) \nonumber\\ &&-\, Q_{d}^{\rm{*}}\,{ +\, \lambda }{\left( {{{\boldsymbol m}} - {{{\boldsymbol m}}^{r}}} \right)^{T}}{{\boldsymbol W}}_{m}^{T}{{{\boldsymbol W}}_{m}}\left( {{{\boldsymbol m}} - {{{\boldsymbol m}}^{r}}} \right),\end{eqnarray} (9)is to be minimized with regard to the set of model parameters contained in the model vector $${{\boldsymbol m}} = {\rm{\ }}{( {{m_1}, \ldots ,{\rm{\ }}{m_{m}}} )^{T}}$$. $${{\boldsymbol d}} = {\rm{\ }}{( {{d_1}, \ldots ,{\rm{\ }}{d_{N}}} )^{T}}$$ is a vector containing N field data points (including RMT and ERT data), and the vector $${{\boldsymbol F}}[ {{\boldsymbol m}} ]$$ contains the corresponding forward responses computed for a given model m. For RMT, the responses are apparent resistivity and phase, and for ERT, the response is apparent resistivity or resistance (for our application, we chose resistance as will be presented later). The superscript T denotes matrix transposition. $${{{\boldsymbol W}}_{d}} = {\rm{\ diag}}{( {\sigma _1^{ - 1}, \ldots ,{\rm{\ \ }}\sigma _{N}^{ - 1}} )^{T}}$$ is a data weighting matrix, where σi is the standard deviation of the observed data. $$\mathbf{W}_m^T\mathbf{W}_m = \alpha_y\mathbf{\boldsymbol\partial}_y^T\mathbf{\boldsymbol\partial}_y + \alpha_z\mathbf{\boldsymbol\partial}_z^T\mathbf{\boldsymbol\partial}_z$$ is the model regularization operator, which controls the simplicity of the inversion model $${{\boldsymbol m}}$$ and contains vertical and horizontal smoothness operators ∂y and ∂z, respectively, with manually adjustable weights αy and αz (Kalscheuer et al.2010). $${{{\boldsymbol m}}^{r}}$$ is a reference model, which is constructed from a priori information. The Lagrange multiplier λ balances data fit and model simplicity. $$Q_{d}^{\rm{*}}$$ is the target data misfit. We use the following definitions of data misfit Qd and root mean square error (RMS) (Kalscheuer et al.2013):   \begin{eqnarray}{Q_{d,sw}}\left[ {{\boldsymbol m}} \right] \!=\! \frac{{{N_d}}}{{\sum\nolimits_{j \!=\! 1}^{{N_{ds}}} {\sum\nolimits_{i \!=\! 1}^{{N_j}} {{{\left( {\frac{1}{{{w_{ji}}}}} \right)}^2}} } }}{\sum\nolimits_{j \!=\! 1}^{{N_{ds}}} {\sum\nolimits_{i \!=\! 1}^{{N_j}} {\left( {\frac{1}{{{w_{ji}}}}\frac{{{d_{ji}} \!-\! {F_{ji}}\left[ {{\boldsymbol m}} \right]}}{{{\sigma _{ji}}}}} \right)} } ^2}\nonumber\\ \end{eqnarray} (10)and   $${\rm{RMS}} = \sqrt {\frac{{{Q_{d,sw}}\left[ {{\boldsymbol m}} \right]}}{{{N_d}}}},$$ (11)where Nd is the total number of data points, Nds is the number of data sets, Nj is the total number of data points in data set j, wji is data set weighting factor and σji is the standard deviation of dji. By linearizing the forward operator in the vicinity of the model of iteration k, $${{{\boldsymbol m}}_{k}}$$, the original problem of minimizing the function $$\Phi ( {{{\boldsymbol m}},{\rm{\lambda }}} )$$ is simplified to minimizing a function $${\Phi ^{{\rm{quad}}}}( {{{\boldsymbol m}},{\rm{\lambda }}} )$$ which is quadratic in $${{{\boldsymbol m}}_{{k} + 1}}$$ (Menke 1989),   \begin{eqnarray} &&{ {\Phi ^{{\rm{quad}}}}\left( {{{{\boldsymbol m}}_{{k} + 1}},{\rm{\lambda }}} \right)}\nonumber\\ &&{\quad = \left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{{{\boldsymbol m}}_{k}}} \right] - {{\bf J}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}_{k}}} \right)} \right)^{T}}\nonumber\\ &&{\qquad \times \,{{\boldsymbol W}}_{d}^{T}{{{\boldsymbol W}}_{d}}\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{{{\boldsymbol m}}_{k}}} \right] - {{\bf J}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}_{k}}} \right)} \right)}\nonumber\\ &&{\qquad - Q_{d}^{\rm{*}}{\rm{\, +\, \lambda }}{\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}^{r}}} \right)^{T}}{{\boldsymbol W}}_{m}^{T}{{{\boldsymbol W}}_{m}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}^{r}}} \right) ,}\end{eqnarray} (12)where $${\boldsymbol{J}} = {\{ {\frac{{\partial {F_{\boldsymbol{i}}}[ {{{{\boldsymbol m}}_k}} ]}}{{\partial {m_j}}}} \}_{m = {m_k}}}$$ is the Jacobian matrix of partial derivatives consisting of N rows and M columns. Thus, i = 1, …, N and j = 1, …, M. In this work, the ERT computation of forward responses and sensitivities in EMILIA has been modified to allow for fields generated by subsurface electrodes. Two types of wavenumber selections for Fourier transformation along the strike direction and inverse Fourier transformation can be chosen in EMILIA (Xu et al.2000; Xiong & Wang 2011). $${{\boldsymbol F}}[ {{{{\boldsymbol m}}_{k}}} ]$$ has been changed to represent either resistances or apparent resistivities, because apparent resistivity is difficult to be defined meaningfully when electrodes are half on land and half under water. The corresponding formulations of sensitivities are presented in Appendix  A. The RMT modeling part in EMILIA did not require any modification for our study. Model error and resolution analysis is required to study model stability and closeness of an estimated model to the true model. The ideal model resolution matrix is an identity matrix of size M × M, which means all model parameters are perfectly resolved (Menke 1989). Generally, the estimated model parameters are weighted averages of the true ones. Thus, the row of the model resolution matrix that pertains to the cell under investigation has non-zero entries also off its diagonal entry (referred to as spread; Menke 1989). Typically, the entries of a row of the model resolution matrix are scaled by the respective cell dimensions yielding a so-called resolving kernel. 3 SYNTHETIC TEST 3.1 ERT forward response: comparison of EMILIA and RES2DMOD In order to evaluate our modification of the ERT part in EMILIA with electrodes located at arbitrary positions, we compared the forward responses computed by EMILIA and RES2DMOD (Loke 2002) for the case of underwater electrodes. The program RES2DMOD (Loke 1999) is a combination of the classic algorithms described by Dey & Morrison (1979) and Silvester & Ferrari (1996). The model used for comparison is originally from Loke (2002), which was slightly modified to achieve a better comparison. It is a two-layer model with the first layer having a resistivity of 50 Ωm (e.g. water) and a thickness of 100 m. The second layer (i.e. the confining half-space) has a resistivity of 100 Ωm. A 10 m × 10 m wide block with a resistivity of 2000 Ωm is located with its top at 5 m depth below the top of the confining half-space. In Fig. 1(a), only the second layer is shown for a good view of the anomaly in the model. Electrodes (red triangles) with a Wenner Alpha array configuration and a spacing of 1 m were placed at the interface, set as 0 m depth, between the first and second layers. The total length of the profile is 50 m. The Wenner Alpha array, which has the strongest signal strength among the common arrays, is normally used for field surveys and usually just referred to as the ‘Wenner’ array (Edwards 1977). Forward responses were computed using the same mesh in the central part of the model with the anomalous block and electrodes with both RES2DMOD and EMILIA (the discretization in the outer parts of the grid towards the boundary nodes is not shown in RES2DMOD). All the simulated data points (408 in total) show relative differences ((ρres2dmod − ρemilia)/ρres2dmod × 100 per cent) between the two codes of less than 1.8 per cent (Fig. 1b). Following the plotting convention for pseudo-sections of the Wenner Alpha array, the horizontal position and depth of each data point are calculated as the midpoint between the potential electrode and 0.519 times the electrode spacing, respectively (Edwards 1977). Figure 1. View largeDownload slide Synthetic model consisting of two conductive layers (interface at z = 0 m, only deeper layer shown) and a 2-D anomaly in the second layer (Fig. 1a) and relative response differences (Fig. 1b) between ERT Wenner forward responses generated by RES2DMOD and EMILIA for the model in (a). Red triangles indicate electrodes. In (b), the horizontal axis represents the middle position of the two potential electrodes and the vertical axis represents 0.519 times a, which is the distance between two abutting electrodes (Edwards 1977). Note that only over limited areas, such as the anomaly and interface between layers 1 and 2 at z = 0 m, relative differences are as high as 1.8 per cent. Figure 1. View largeDownload slide Synthetic model consisting of two conductive layers (interface at z = 0 m, only deeper layer shown) and a 2-D anomaly in the second layer (Fig. 1a) and relative response differences (Fig. 1b) between ERT Wenner forward responses generated by RES2DMOD and EMILIA for the model in (a). Red triangles indicate electrodes. In (b), the horizontal axis represents the middle position of the two potential electrodes and the vertical axis represents 0.519 times a, which is the distance between two abutting electrodes (Edwards 1977). Note that only over limited areas, such as the anomaly and interface between layers 1 and 2 at z = 0 m, relative differences are as high as 1.8 per cent. 3.2 Inversion: synthetic example for single and joint inversions Before presenting the use of the modified implementation for the field case study, a synthetic example is shown. A test model (Fig. 2a) was abstracted from the single inversion result of the ERT field data (shown later). In the model, very conductive sea water with a resistivity of 1.38 Ωm exists at shallow depth in the central part (160–600 m) of the profile overlying bedrock with a resistivity of 1000 Ωm. A fracture zone (100 Ωm) that extends from 600 to 650 m along the profile and dips to northwest at an angle of about 50° from the horizontal is contained in the model (since the code only can use rectangular mesh, the fracture zone has step-like boundaries). The positions of RMT and ERT stations are marked by red triangles and black dots, respectively (Figs 2b and c). On the lake, the RMT station spacing is around 10 m and, on land, it varies from 20 to 40 m (yielding 52 stations in total). However, the frequency range is fixed from 14 to 250 kHz (nine frequencies equally distributed in log space). ERT data were generated for a gradient array type. The spacing between two adjacent electrodes is 5 m. The total length of the ERT cable is 790 m (as for the field example, data from 156 electrodes were used and data from one electrode were discarded owing to coupling problems or instrumental malfunctioning). EMILIA was used to model the RMT responses on the surface including water and land data and the ERT responses on land and at the bottom of the water body. Then, we contaminated the synthetic data with Gaussian noise of 5 per cent on ERT resistance, 10 per cent on RMT apparent resistivity and 2.29° on RMT phase data, which are on levels comparable to those observed in our field data. Note that the noise levels selected for the RMT apparent resistivities and phases correspond to a 5 per cent noise level on the impedance tensor elements. This corresponds to 5 per cent on ERT resistance data, if we assume that the errors in the RMT measurements are mainly in the electrical field. An experiment was carried out to compare the differences between individual and joint inversions with the synthetic data sets. Figure 2. View largeDownload slide (a) Test model deduced from single inversion of ERT field data (see below). (b) Single inversion model of RMT TE-mode data. Red triangles represent RMT stations at the surface and white dashed line shows the approximate exploration depth of RMT signals. (c) Single inversion model of ERT data. Black dots represent ERT electrode positions and white dashed line shows the approximate exploration depth of ERT. (d) Joint inversion model for RMT TE-mode and ERT data. White dashed line shows the approximate joint exploration depth of RMT and ERT signals. (e)Table showing RMS misfits. The joint inversion model can fit both data sets equally well. Figure 2. View largeDownload slide (a) Test model deduced from single inversion of ERT field data (see below). (b) Single inversion model of RMT TE-mode data. Red triangles represent RMT stations at the surface and white dashed line shows the approximate exploration depth of RMT signals. (c) Single inversion model of ERT data. Black dots represent ERT electrode positions and white dashed line shows the approximate exploration depth of ERT. (d) Joint inversion model for RMT TE-mode and ERT data. White dashed line shows the approximate joint exploration depth of RMT and ERT signals. (e)Table showing RMS misfits. The joint inversion model can fit both data sets equally well. Only the joint inversions of RMT TE-mode and ERT data are shown in this paper. For the RMT TM mode and the ERT method, the directions of current flow in a 2-D model are relatively similar with current flowing exclusively and predominantly, respectively, in the plane of the profile (i.e. perpendicular to the strike direction in both horizontal and vertical directions). Thus, including RMT TM-mode and ERT data in a joint inversion may lead to inclusion of relatively similar model constraints. Hence, joint inversion of RMT TE-mode and ERT data, where the current flows are perpendicular and predominantly parallel (owing to the use of point sources, there is an additional component of ERT current flow along the strike direction) to the plane of the profile, respectively, leads to more complementary model constraints than joint inversion of RMT TM-mode and ERT data. However, RMT TM-mode and determinant data results are shown in Appendix  B together with a discussion of the preferred use of the TE-mode data in the joint inversions. The single inversion model of the RMT TE-mode data is shown in Fig. 2(b). The shape of the water body is resolved. However, the fracture zone is not clearly resolved. The individual inversion of ERT data gives a similar result (Fig. 2c). White dashed lines represent the maximal exploration depths of the RMT and ERT data sets in the models estimated using Spies’ (1989) method and Dahlin and Zhou's (2006) method, respectively. The maximal RMT exploration depth was estimated using the lowest signal frequency of 14 kHz and the vertical resistivity section underneath each station. However, rather than using 1.5 skin depths as proposed by Spies (1989) to estimate the depth to which structure can be detected, we employed one skin depth to retrieve a conservative estimate for the depth to which the model is well constrained by the data. The joint inversion of RMT TE-mode and ERT data generated a model very similar to the true model (Fig. 2d). Both water and fracture zones are more accurately reconstructed than in any individual inversions. Especially, the dipping angle of the fracture zone is clearly evident. The contrast between the dipping conductor at 150 m distance and its surroundings is not particularly well pronounced. We would not interpret it as a fracture zone, and this feature may mainly be caused by the noise we added to our synthetic data. 4 CASE STUDY In the frame of the TRUST (TRansparent Underground STructure; TRUST 2013) project, boat-towed RMT and lake-floor ERT data were collected. The main aim of the project is to develop and implement methods and techniques to obtain more accurate models for planning, designing and constructing urban underground transportation systems, such as tunnels, subways and trams (Bastani et al.2015; Brodic et al.2015; Malehmir et al.2015; Brodic et al.2017; Mehta et al.2017). Äspö was chosen as one of several common sites for testing the effectiveness of the combined use of geophysical, geochemical and hydrological methods. Boat-towed RMT (Bastani et al.2015) and lake-floor ERT surveys (Dahlin et al.2014), among other methods, were carried out in May 2015 with the objective of imaging geological structures below a lake individually and in combination. 4.1 Geological setting Äspö is located approximately 30 km north of the Oskarshamn archipelago near the shoreline of the Baltic Sea in southeastern Sweden (Fig. 3). Granitic rocks and diverse types of fractures or fracture zones dominate the site (Cosma et al.2001). Most fracture zones at Äspö are a result of reactivation of older structures and appear mainly brittle. The style of such fracture zone tends to depend on the nature of any older structure being reactivated (Stanfors et al.1999). The lake where we performed our measurements is connected to the Baltic Sea through a narrow water channel. In the northernmost part of the site, the Äspö HRL operated by SKB (Swedish Nuclear Fuel and Waste Management Company) is located. A nuclear power plant is about 500 m away from the surveyed line. An underground tunnel (partly shown as an NNW trending blue solid line in Fig. 3b) plunging at about 14 per cent over a length of about 1500 m connects the surface to the HRL as well as various smaller tunnels, ramps and one main shaft (Almén & Stenberg 2005). Figure 3. View largeDownload slide (a) Areal photo (Google map) and (b) geological map of the case study area (courtesy of SGU). RMT profiles are marked by different symbols with different colours. The RMT profile used for the joint inversion is marked by red circles on the lake and pink triangles on the land. ERT data were measured only for one profile, marked by green dots. Positions of fracture zones, which are documented by SKB, are shown in (b). However, they are projected to the ground surface based on limited tunnel and borehole observations combined with lost low-resolution geophysical data. The bedrock is mainly granitic. The general topography of the land part is approximately flat. Figure 3. View largeDownload slide (a) Areal photo (Google map) and (b) geological map of the case study area (courtesy of SGU). RMT profiles are marked by different symbols with different colours. The RMT profile used for the joint inversion is marked by red circles on the lake and pink triangles on the land. ERT data were measured only for one profile, marked by green dots. Positions of fracture zones, which are documented by SKB, are shown in (b). However, they are projected to the ground surface based on limited tunnel and borehole observations combined with lost low-resolution geophysical data. The bedrock is mainly granitic. The general topography of the land part is approximately flat. Based on the information obtained during tunnel construction, an NE-SW running fracture system (known as NE-1) exists below the northern end of the RMT lake profile (Fig. 3b). It can be regarded as a southern boundary to the laboratory and plays an important role for the geometries of the fractures and deformation zones that formed at the Äspö HRL site (Berglund et al.2003). At about 1300 m along the access tunnel, the NE-1 fracture system is intersected at approximately 180 m below ground surface (Almén & Stenberg 2005). Three main subsets comprise the fracture system and are about 60 m wide in total (Berglund et al.2003; Rhén et al.1997). The two southernmost subsets, trending NE and dipping to NW, can be described as highly fractured and hydraulically conductive zones (Stanfors et al.1999; Makurat et al.2006). The northern subset, approximately 28 m wide, is the most intensely fractured part of NE-1 and produced significant inflow of water to the tunnel during the construction phase (Stanfors et al.1999; Makurat et al.2006; Almén & Stenberg 2005). Diorite, fine-grained granite and greenstone are the main host rocks of the zone (Stanfors et al.1999; Berglund et al.2003). The central part of NE-1 is strongly altered to clay. Since NE-1 is strongly water bearing, various other fractures around the zone with various orientations (Fig. 3b) are also water bearing (Berglund et al.2003). The water is a mixture of non-saline and brackish sea water (Wikberg et al.1991). The dip of the NE-1 fracture zone in the tunnel intersection and adjusted by cored boreholes is about 65° with uncertainty, because information gaps exist between the tunnel and the boreholes (Berglund et al.2003). Further away from the tunnel and boreholes, the complex geological characteristics of the zone are only poorly determined. Furthermore, the NE-1 fracture zone is not exposed at the surface. Hence, geophysical information becomes especially important for resolving its extent and geometry (Berglund et al.2003). Three fracture zones namely EW-7, NE-4 and NE-3 in the southern half of the RMT profile below the lake (Fig. 3b) were indicated by refraction seismic (not publically available) and borehole data (Wikberg et al.1991; Stanfors et al.1999; Rhén et al.1997). The NE-4 fracture zone trends NE and consists of two continuous subsets with an approximate width of 40 m. In the access tunnel, the NE-3 fracture zone is 49 m wide. It dips steeply towards NNW (Stanfors et al.1999). Geological information suggests the existence of the fracture zone EW-5 (Fig. 3b, Wikberg et al.1991); however, it cannot be observed in the tunnel. In our study, new geophysical data along with bathymetry data and water resistivity measurements are used in 2D inversions to delineate the shape of the NE-1 fracture zone and others that are inferred such as EW-5. 4.2 Data acquisition and quality The instrument Enviro-MT (Bastani 2001) developed at Uppsala University (UU) was used for RMT data acquisition. This system was upgraded to a boat-towed data acquisition system in 2015 thanks to the collaboration between UU and the Geological Survey of Sweden (SGU, Bastani et al.2015). In our field measurements, the RMT station spacing was not fixed. On the lake the spacing was around 10 m, but on land the spacing varied from 20 to 40 m depending on where sufficient space was available to place the RMT system (Fig. 3a) resulting in 52 RMT stations in total. The signal frequency range of the received Very Low Frequency and Low Frequency radio transmitters is from 14 to 250 kHz, which as a result of further processing following Bastani & Pedersen (2001) gave nine frequencies equally distributed in log space. Compared to the land part, the RMT data acquisition on the lake was very efficient. Within two and a half hours, all the RMT stations on the lake were surveyed (Fig. 3a). However, two days were used to acquire the land RMT data (Fig. 3a). RMT stations surveyed and used for joint inversion are marked by red circles on the lake and pink triangles on land. The instrument Terrameter LS, designed by ABEM, was used to carry out ERT data acquisition using a multiple gradient array (Dahlin & Zhou 2006). The spacing between two adjacent electrodes was 5 m and the ERT electrode spread was 635 m long. The survey line was extended by roll-along to a total length of 790 m (after editing, data from 156 electrode positions out of a total of 157 electrode positions could be used). However, since the multielectrode cable was placed along the lake floor combined with some parts on land, the horizontal distance is 783.5 m along the profile. Along the lake floor, the electrode outlets of the cable were placed in the sediments (the electrode cable sank into the sediments when it was unrolled from a boat, because its density is higher than that of water). Since this provided sufficiently low coupling resistance, connecting electrodes to the cable was rendered unnecessary. The cable outlets on land were connected to plate electrodes, which were buried where glacial sediments were present and attached on the granite bedrock with a conductive gel on outcrops. More details of the ERT data acquisition can be found in Ronczka et al. (2017). ERT stations used for joint inversion are marked by green dots in Fig. 3(a). The field data sets are shown in Fig. 4. The RMT TE-mode apparent resistivity and phase data collected on water suggest that the sea water (characterized by the signals with frequencies >100 kHz) has even higher resistivity than the underlying sediments (characterized by the 50–100 kHz signals; Figs 4a and b), which together with the sea water strongly reduced the penetration depth of the RMT signals (eq. 3). The total thickness of water and sediments is small (<3.5 m) along 50 per cent of the part of the profile covered by water, for example, at 380–600 m distance along the profile (Fig. 4d), and around 20 per cent of the RMT stations are on land. Therefore, the RMT data can provide useful information for modeling the underlying resistive bedrock and the lake sediments. Our vertical magnetic transfer function data are not of sufficiently high quality to be used in the interpretation (convergence is bad when inverting magnetic transfer function data). Since the ERT cable was placed on the lake floor, the ERT penetration depth in the underlying resistive basement is not affected strongly by the overlying conductive water. The ERT data in Fig. 4(c) show low-resistance values in the centre of the line and high-resistance features at both ends of the line corresponding to outcropping bedrock. Before the data sets were inverted, noisy data were removed. Owing to the diffusive and integrative natures of EM and electrical fields, the criterion for removing data points was deviation from smooth response variation. For the inversion, a 10 per cent relative error floor was used for the RMT apparent resistivity data (30 per cent relative error floor was used for the land RMT apparent resistivity data to avoid influence from static shift), an absolute error floor of 2.29° was used for RMT phase data, and a 5 per cent relative error floor was used for the ERT data. Figure 4. View largeDownload slide Field measurements along the coincident RMT and ERT profiles: (a) RMT TE-mode apparent resistivity data, (b) RMT TE-mode phase data, (c) ERT resistance data and (d) elevation of lake floor. The lake-floor sediments have lower resistivity than water which reflects in RMT apparent resistivity (panel a) at frequencies of 50–100 kHz and at 300–500 m distances along the profile being lower than at the highest frequencies, which penetrate only into water. Figure 4. View largeDownload slide Field measurements along the coincident RMT and ERT profiles: (a) RMT TE-mode apparent resistivity data, (b) RMT TE-mode phase data, (c) ERT resistance data and (d) elevation of lake floor. The lake-floor sediments have lower resistivity than water which reflects in RMT apparent resistivity (panel a) at frequencies of 50–100 kHz and at 300–500 m distances along the profile being lower than at the highest frequencies, which penetrate only into water. 4.3 RMT strike analysis Based on the relevant geological information and previous studies of the site, the RMT data acquisition parameters (profile orientation and station spacing) were carefully chosen to resolve possible fracture zones (Fig. 3). However, dimensionality, distortion and strike analysis of the RMT data needs to be done to decompose the RMT field data into TE- and TM-mode data. A number of publications discussed the strike rules for the magnetotelluric impedance tensor in details (e.g. Zhang et al.1987; Jones & Groom 1993). In our study, the method proposed by Zhang et al. (1987) was used for dimensionality, distortion and strike analysis by calculating distortion parameters and strike directions from the impedance tensor and evaluating the fit of the distortion model to the field data. When applying this method to the RMT data, independent strike angles are obtained for each station and frequency. For our data, only one preferred direction of strike shows in a cumulative rose diagram (Fig. 5a), which is 83°E. Since we oriented the sensors of both the land and boat-towed RMT stations parallel and perpendicular to the profile direction, this strike direction is relative to the profile direction, and it corresponds to 75° East of geographic North. Swift (1967) skews of all of our RMT data (Fig. 5b, except for the impedance tensor of one station at the lowest signal frequency) are lower than 0.2 and most of them are lower than 0.1. Therefore, our RMT data can be considered to represent a 1-D/2-D structure and 3-D effects should not become problematic in 2-D inversions and in later interpretation of the 2-D resistivity model. According to Fig. 3(b), the RMT profile used for joint inversion is almost perpendicular to the speculated fracture zone orientation. Since the field measurements were well designed using the available fracture zone information and the joint inversion of RMT and ERT data requires coincident profile directions, it turned out that there was no need for rotation of the reference coordinate system to the strike direction. However, projection of the first five station locations onto a modeling profile perpendicular to the strike direction at the southeastern end was needed, because these stations are 35–50 m away from the profile (Fig. 3a). Figure 5. View largeDownload slide (a) Strike directions calculated with the method proposed by Zhang et al. (1987). (b) Swift skews. Both of them show the structure underneath the research area is 2-D with a preferred strike direction of 83° East with regard to the profile direction used in the field (corresponding to 75° East of geographic North). Figure 5. View largeDownload slide (a) Strike directions calculated with the method proposed by Zhang et al. (1987). (b) Swift skews. Both of them show the structure underneath the research area is 2-D with a preferred strike direction of 83° East with regard to the profile direction used in the field (corresponding to 75° East of geographic North). 4.4 Single and joint inversion results We carried out the inverse modeling in two steps. In the first step, we performed an Occam inversion using regular smoothness constraints (Constable et al.1987; Menke 1989) with a variable Lagrange multiplier and a 100 Ωm half-space initial model. In the second step, the final inversion model from step one was then used as a new initial model in an Occam inversion with additional Marquardt–Levenberg damping. Here, the Langrage multiplier for the smoothness constraint is kept fixed at the value that gave the lowest RMS in step one, while an optimal damping factor (values varied from 100.0 to 102.8 in our inversions) is determined in each iteration using a line search. Note, in a few cases using the average resistivity of the final inversion model from step one as a new initial model for step two resulted in a better data fit. Consequently, we considered this model for further considerations. The models from individual and joint inversions of the field data are shown in Fig. 6. Both individual inversions show a conductive water body, more conductive sediments saturated with saline water below the water body (ERT electrodes positions mark the water bottom in Fig. 6b) and a high-resistivity feature below the lake floor in the central part of the profile (at approximately 400 m along the profile) most possibly correlating with the small island east of the profile (Fig. 3). Besides, two low-resistivity zones (possibly fracture zones) appear at distances of 200–350 and 550–700 m along the profile in both single inversion models (Figs 6a and b). Comparing the two models (Figs 6a and b), the RMT data provide more information of the possible fracture zone at 550–600 m distance along the profile than the ERT data (white dashed lines indicate exploration depth in both models). However, there is a strong similarity between the conductive structures in both models. The joint inversion model (Fig. 6c) shows two relatively conductive zones at 550–700 m distance along the profile. Both of them are above the exploration depth (the white dashed line, which is the maximal exploration depth of the joint RMT and ERT data sets). The same consideration of exploration depth is suitable for the conductive zone at 200–350 m distance along the profile. The RMS misfits (in panel d of Fig. 6) of the joint and single inversion models are comparable; however, the joint inversion model can fit two different geophysical data sets simultaneously and more detailed structures are introduced in the model. This is because the null space (i.e. the part of the model space which is not determined by the data, e.g. Menke 1989) corresponding to the RMT TE data can partly be resolved by the ERT data and vice versa. This is also well proven by another experiment as follows. The models from single inversions of RMT and ERT data were used as initial models for single inversions of the ERT and RMT data, respectively, providing initial RMS misfits of 9 in either case. Considering the final RMS of 2.9 in the joint inversion (Fig. 6) implies that a single inversion model from one data set does not explain the other data set well, because each data set has its own limitations. In such a situation, joint inversion plays a key role to reduce the null space due to the differences in data coverage, sensitivity and noise influences. Figure 6. View largeDownload slide Inversion of field data: (a) Inversion model computed from RMT TE-mode data. (b) Inversion model computed from ERT data. (c) Joint inversion model for RMT TE-mode and ERT data. (d) Table showing RMS. Separate RMS values from joint inversion are slightly higher than those of the single inversions. However, the model fits both data sets with acceptable RMS. Figure 6. View largeDownload slide Inversion of field data: (a) Inversion model computed from RMT TE-mode data. (b) Inversion model computed from ERT data. (c) Joint inversion model for RMT TE-mode and ERT data. (d) Table showing RMS. Separate RMS values from joint inversion are slightly higher than those of the single inversions. However, the model fits both data sets with acceptable RMS. Although TE-mode data were chosen for this study, TM, TE and TM, and determinant data were also inverted for 2-D models. In Appendix  B, single inversion models of TM-mode and determinant impedances and joint inversions combined with ERT data are presented (see also Fig. B3). 4.5 Weighting of data sets and smoothness constraints ERT and RMT data sets usually have different data point densities and different sensitivities to model structures and are affected differently by 3-D anomalies implying the need of different weights on the ERT and RMT data sets (Kalscheuer et al.2013; Sudha et al. 2014). Since the ERT data set has more data points and better data coverage than the RMT data set, more weight on the ERT data set only results in the joint inversion model being close to the ERT single inversion model (Fig. 7a). The RMS values of the individual data sets also prove that the joint inversion is dominated by the ERT data (see the caption of Fig. 7). Moreover, higher weight on the RMT data set leads to inclusion of more information from the RMT data, which is superior to resolve conductors (in our case the possible fracture zone). Joint inversion with two times more weight on the RMT data set shows two conductive zones at 550–700 m distance along the profile (Fig. 7b). Those two conductive zones are consistent with the ones in the inversion (Fig. 6c) without any weighting. Other weights (1.5:1, 3:1 and 4:1) were also tried without giving better results than the one shown. Figure 7. View largeDownload slide (a) Joint inversion model using weight on ERT data two times higher than on RMT data. Apparently, this inversion is dominated by ERT data set. (b) Joint inversion model using weight on RMT data two times higher than on ERT data. Since the ERT data set has a higher number of data points, higher weight on the RMT data leads the joint inversion to adequately fit both data sets. (c) Joint inversion model when weight on horizontal model smoothness is twice the weight on vertical smoothness. (d) Joint inversion model constrained with bathymetry data and water resistivity measurement. (e) Zoom-in figure of part of the model in (d) marked by a rectangle. (f) Table showing RMS. All symbols in the subfigures are the same as in Fig. 6. All joint inversions are consistent with regard to the conductive zones along the profile except for the ERT dominated one. Figure 7. View largeDownload slide (a) Joint inversion model using weight on ERT data two times higher than on RMT data. Apparently, this inversion is dominated by ERT data set. (b) Joint inversion model using weight on RMT data two times higher than on ERT data. Since the ERT data set has a higher number of data points, higher weight on the RMT data leads the joint inversion to adequately fit both data sets. (c) Joint inversion model when weight on horizontal model smoothness is twice the weight on vertical smoothness. (d) Joint inversion model constrained with bathymetry data and water resistivity measurement. (e) Zoom-in figure of part of the model in (d) marked by a rectangle. (f) Table showing RMS. All symbols in the subfigures are the same as in Fig. 6. All joint inversions are consistent with regard to the conductive zones along the profile except for the ERT dominated one. Effects of different weights on the horizontal and vertical smoothing of the model (Kalscheuer et al.2010) were also evaluated using weights of one on both data sets. More weight on horizontal smoothing was used to generate structures that are extended predominantly in the horizontal direction. Twice more weight on horizontal smoothing than on the vertical smoothing did not smear out any of the vertical conductive zones at 550–700 m distance along the profile (Fig. 7c) suggesting that these conductive zones are required by the data. In this test with increased weight on horizontal smoothing, using the data set weights that were previously found to be optimal (i.e. two times more weight on the RMT data than on the ERT data) resulted in increased RMS misfits (total RMS = 3.61, RMT TE-mode data RMS = 3.05 and ERT data RMS = 4.32). This seems to be related to the fact that the RMT TE-mode data and ERT data have different sensitivities to horizontal and vertical resistivity contrasts thus necessitating adjustment of data set weights for changes in model smoothness constraints and vice versa. 4.6 Joint inversion using an a priori model During the RMT lake measurement, bathymetry data were also acquired. The water resistivity was measured at three different depths (Ronczka et al.2017). A three-layer model based on the water resistivity measurements (1.48 Ωm at 0.2 m depth, 1.37 Ωm at 3.8 m depth and 1.28 Ωm at 4.8 m depth.) can be used to represent the resistivity of the lake water. This information was used as additional constraints in the joint inversion by constructing an a priori model (mr in eq. 12). The a priori model had two structures, the water body with three layers and the rest of the model consisting of bedrock with a resistivity of 225 Ωm. The geometry of the low boundary of the water body in the a priori model was controlled by the bathymetry data. The bedrock resistivity in the a priori model was the average apparent resistivity value of presumed bedrock structure in the joint inversion models in Figs 6(c) and 7(a)–(c). Since the water resistivity was known at three different depths, it was fixed in the hope that this would constrain the remaining model parameters to find representative values. However, below the lake floor, more conductive sediments were observed in the inversion models without bathymetric constraints. Hence, in implementing the bathymetric constraints, the smoothness constraints along the lake floor were removed to guarantee the free change of the model parameters across the lake floor. By this way, the water resistivity and the lake floor depth were introduced in the inversion without influencing the parts of the model lacking the a priori information. The two-step inversion strategy presented in Section 4.4 was used. The initial model in step one is the a priori model built from the available information, and the a priori model and smoothness constraints constructed according to water resistivity, average apparent resistivity and bathymetry were used in both inversion steps. The joint inversion model of the RMT TE-mode and ERT data constrained by bathymetry and measurements of water resistivity (Fig. 7d) is similar to the results without these constraints. In this model, the shape of the lake-floor sediments is more focused (Fig. 7e). Also, the bedrock surface is better defined (as compared to the other joint inversions) and the resistivity values of bedrock are more representative of crystalline rocks (mostly above 10 000 Ωm). The conductive zones at 550–700 m distance along the profile are at least as well resolved as in the other joint inversions. In this inversion, the weights of the RMT TE-mode and ERT data sets and of the horizontal and vertical smoothness constraints were set to one. In alternative inversions using weights deviating from one (results not shown), the RMS misfits were higher. 4.7 Evaluation of inversion models We evaluate the joint inversion results by analysing (1) the exploration depth (already shown above), (2) the data misfit for each data point and (3) model resolution and error estimates. Data misfits ((dfield − demilia)/σfield, where σfield represents data uncertainties) are displayed for all the inversion models in Figs 8 and 9. In general, the RMT phase data have better fits than the RMT apparent resistivity data (Fig. 8). Compared with the RMT data misfit values (Fig. 8), the ERT data misfits are lower (Fig. 9). Importantly, joint inversion models without bathymetric and water resistivity constraints show data misfits comparable to those of the single inversion results (Figs 8 and 9). In particular for the RMT data, the joint inversion with bathymetric and water resistivity constraints shows higher RMS than the single inversions. However, a careful look at the data differences suggests that the higher frequency data corresponding to the water layer apparently show higher misfits (Fig. 8f). This is likely because the assumption of a three-layer water resistivity model is not accurate enough. A one-layer model was tested using an average of the water resistivities and a complex model was also tested using a linear decrease of the water resistivity with depth, but both inversions showed slightly degraded fits compared to the three-layer model. The conductive zones at 500–700 m distance are better resolved in the joint inversion model computed with constraints from bathymetry and water resistivity (Fig. 7d) yielding small data differences comparable to the single inversion results. This feature can be observed in Fig. 8(f) and Table 1, especially in the fit to the RMT apparent resistivity. Since the ERT data were collected on the lake sediments, the likely inaccurate assumption of water resistivity seems to have little influence and a careful look at the misfit reveals the improvement of data misfit (the elliptic zone marked in Fig. 9f as compared to Figs 9a–e). Thus, the joint inversion result constrained by bathymetry and water resistivity in Fig. 7(d) was chosen for interpretation. In general, the joint inversion improves the details of the inverted model and results in fitting both data sets with RMS values similar to those of single inversions. Figure 8. View largeDownload slide Normalised data differences between field and inverse modeled RMT data from (a) RMT TE data single inversion; (b) joint inversion without any weighting and a priori information; (c) joint inversion when weight on ERT data is two times higher than on RMT data; (d) joint inversion when weight on RMT data is two times higher than on ERT data; (e) joint inversion when weight in horizontal direction on the model is two times higher than in vertical direction and (f) joint inversion when constrained with bathymetry data and water resistivity measurements. Most joint inversions have RMS misfits comparable to the single inversion. In general, the phase is better fitted than the apparent resistivity. Crosses show the data points removed before the inversions were run (around 13 per cent of apparent resistivity data and 16 per cent of phase data). In the constrained joint inversion (f), a relatively large misfit is observed, because describing the water column with three layers is not accurate enough to represent the spatial variability of water resistivity. However, the part below the water zone gives comparable or even better [apparent resistivity in (f)] misfit compared to single inversions. The relevant RMS values are shown in Table 1. Figure 8. View largeDownload slide Normalised data differences between field and inverse modeled RMT data from (a) RMT TE data single inversion; (b) joint inversion without any weighting and a priori information; (c) joint inversion when weight on ERT data is two times higher than on RMT data; (d) joint inversion when weight on RMT data is two times higher than on ERT data; (e) joint inversion when weight in horizontal direction on the model is two times higher than in vertical direction and (f) joint inversion when constrained with bathymetry data and water resistivity measurements. Most joint inversions have RMS misfits comparable to the single inversion. In general, the phase is better fitted than the apparent resistivity. Crosses show the data points removed before the inversions were run (around 13 per cent of apparent resistivity data and 16 per cent of phase data). In the constrained joint inversion (f), a relatively large misfit is observed, because describing the water column with three layers is not accurate enough to represent the spatial variability of water resistivity. However, the part below the water zone gives comparable or even better [apparent resistivity in (f)] misfit compared to single inversions. The relevant RMS values are shown in Table 1. Figure 9. View largeDownload slide Normalised data differences between field and modeled ERT data for the same [except for the single inversion in (a)] inversion models as considered in Fig. 8. Joint inversions have RMS misfits comparable to the single inversion. The data marked by an ellipse at 600 m along the profile cause high misfit. Apparently, the joint inversion constrained by bathymetry and measurements of water resistivity shows the best fit over this area. Figure 9. View largeDownload slide Normalised data differences between field and modeled ERT data for the same [except for the single inversion in (a)] inversion models as considered in Fig. 8. Joint inversions have RMS misfits comparable to the single inversion. The data marked by an ellipse at 600 m along the profile cause high misfit. Apparently, the joint inversion constrained by bathymetry and measurements of water resistivity shows the best fit over this area. Table 1. Averages of data differences shown in Fig. 8. However, the three highest frequencies mainly relating to saline water are excluded in order to have a better comparison among the different types of inversion (cf. main text). The constrained joint inversion (Fig. 8f) has obviously the smallest data differences compared with the single inversion (Fig. 8a). Data differences  Fig. 8(a)  Fig. 8(b)  Fig. 8(c)  Fig. 8(d)  Fig. 8(e)  Fig. 8(f)  Average ρ difference  2.08  2.59  2.78  2.43  2.83  1.93  Average φ difference  0.79  1.16  1.91  0.94  1.26  1.24  Data differences  Fig. 8(a)  Fig. 8(b)  Fig. 8(c)  Fig. 8(d)  Fig. 8(e)  Fig. 8(f)  Average ρ difference  2.08  2.59  2.78  2.43  2.83  1.93  Average φ difference  0.79  1.16  1.91  0.94  1.26  1.24  View Large A model error and resolution analysis was then done for the joint inversion model constrained by bathymetry and water resistivity in Fig. 7(d). A model resolution analysis determines to what extent the true model maps into an inversion model. The smaller the spread of non-zero entries in the row of the resolution matrix around the diagonal entry is and the higher the diagonal entry is, the better is the model resolution provided by the data (Kalscheuer & Pedersen 2007; Kalscheuer et al.2010, 2013). Due to the limited space, the resolving kernels of only a few cells of the joint inversion model constrained by bathymetry and water resistivities are shown in Fig. 10 and the relevant information of the selected cells as well as the linearized model errors f are listed in Table 2. Underneath the land parts of the profile, the model resistivities are partly resolved to a depth of about 20 m (e.g. cell a in Fig. 10a). However, the main lobe of the resolving kernel is already vertically offset from the investigated cell. The sediments below the lake are also rather well resolved (cell b in Fig. 10b), but below the sediments the corresponding resolving kernels are not focused (cells c, e and d in Figs 10c and d) because the conductive water and sediments together limit the penetration depth considerably. Cells representative of the speculated fracture zones are partly resolved (cells f and g in Figs 10a and b). However, the corresponding resolving kernels show only positive side lobes on the cells under investigation while the main lobes are located at shallow depth (<10 m). Possible reasons for the limited resolution of the speculated fracture zones are noise in the data (the RMT and ERT data are noisy in those parts, see Figs 8 and 9), the general experimental setup (one of our main targets is located close to the NW end of the profile) and fundamental physical equivalences (non-uniqueness of the ill-posed inverse problem). Figure 10. View largeDownload slide Resolving kernels for seven model parameters (a–g) of the joint inversion model constrained by bathymetry and water resistivity measurements in Fig. 7(d). The considered cells are marked by white diamonds in the figure panels. The red crosses depict the centres of resolution and the horizontal and vertical resolution lengths (Kalscheuer & Pedersen 2007). The isolines are the logarithmic resistivity values of the model in Fig. 7(d). Two cells are shown in one subfigure to take less space. Model parts with low resolution show a shift between the investigated cells (diamonds) and the main lobes of the resolving kernels (cell d as an example). Figure 10. View largeDownload slide Resolving kernels for seven model parameters (a–g) of the joint inversion model constrained by bathymetry and water resistivity measurements in Fig. 7(d). The considered cells are marked by white diamonds in the figure panels. The red crosses depict the centres of resolution and the horizontal and vertical resolution lengths (Kalscheuer & Pedersen 2007). The isolines are the logarithmic resistivity values of the model in Fig. 7(d). Two cells are shown in one subfigure to take less space. Model parts with low resolution show a shift between the investigated cells (diamonds) and the main lobes of the resolving kernels (cell d as an example). Table 2. Positions (yc, zc) and extents (Δy, Δz) of the cells chosen for model error and resolution analyses as well as linearized error factors f of the resistivities of these cells. Resolving kernels are depicted in Fig. 10. Cell lable  yc  zc  Δy  Δz  f  a  100.2  13.2  4.9  2.1  1.09  b  199.6  18.5  5.0  3.4  1.10  c  299.4  22.3  5.0  4.3  1.10  d  449.1  27.2  5.0  5.4  1.10  e  499.1  22.3  5.0  4.3  1.10  f  548.7  18.5  5.0  3.4  1.13  g  597.7  13.2  5.0  2.1  1.08  Cell lable  yc  zc  Δy  Δz  f  a  100.2  13.2  4.9  2.1  1.09  b  199.6  18.5  5.0  3.4  1.10  c  299.4  22.3  5.0  4.3  1.10  d  449.1  27.2  5.0  5.4  1.10  e  499.1  22.3  5.0  4.3  1.10  f  548.7  18.5  5.0  3.4  1.13  g  597.7  13.2  5.0  2.1  1.08  View Large 4.8 Geological interpretation The geological interpretation is summarized in Fig. 11(a). Three distinct fracture zones exist at distances of 200–300 m along the profile based on previous studies (Wikberg et al.1991; Stanfors et al.1999). These fracture zones are not well resolved in our models (EW-7 is the only one detected, Fig. 11a) because the overlying conductive sea water and sediments limit the resolution in this part of the profile. Following Wikberg et al. (1991) and Stanfors et al. (1999), the fracture zones NE-3, NE-4 and EW-7 should exist in this part of the profile. In the resistivity models shown, the sediments above the fracture zones are more conductive than the sea water (see also Ronczka et al.2017). At 500–700 m distance along the profile, the model suggests two distinct conductive zones (possible fracture zones). The one from 500–580 m distance is likely related to EW-5 (Wikberg et al.1991), which was previously a poorly constrained fracture zone. However, more information is needed to conclusively confirm the existence of EW-5. The other conductive zone at 600–700 m distance along the profile is most probably related to the NE-1 fracture system. The 65° dip of NE-1 suggested by Berglund et al. (2003) matches closely the models obtained by the joint inversion (Fig. 11a). The position of NE-1 in the inversion models matches the surface projection of NE-1 observed in the tunnel (Figs 3b and 11a). However, the projection could be slightly shifted compared with the real position, because it is projected from the tunnel and limited boreholes to the ground surface and no relevant observation is available on the ground surface. Figure 11. View largeDownload slide (a) Interpretation of joint inversion model in Figure 7(d) with possible positions of fracture zones EW-7, NE-4, NE-3, EW-5 and NE-1. Black thick arrows on top of the model section mark the positions of the fracture zones along the profile based on the previous studies shown in Fig. 3. Black lines at depth mark interpreted boundaries of fracture zones as inferred from the model in panel (a). (b) Test model generated with main structures from joint inversion model constrained by bathymetry and water resistivity measurements at 500–700 m distance and from previous studies at 200–350 m distance (EW-7 could be observed in our model). (c) Single inversion model of synthetic RMT TE-mode data. (d) Single inversion result of synthetic ERT data. (e) Joint inversion model of synthetic RMT TE-mode and ERT data. (f) Table showing RMS. The joint inversion model fits both data sets reasonably well and shows the highest level of detail among all the inversions. Figure 11. View largeDownload slide (a) Interpretation of joint inversion model in Figure 7(d) with possible positions of fracture zones EW-7, NE-4, NE-3, EW-5 and NE-1. Black thick arrows on top of the model section mark the positions of the fracture zones along the profile based on the previous studies shown in Fig. 3. Black lines at depth mark interpreted boundaries of fracture zones as inferred from the model in panel (a). (b) Test model generated with main structures from joint inversion model constrained by bathymetry and water resistivity measurements at 500–700 m distance and from previous studies at 200–350 m distance (EW-7 could be observed in our model). (c) Single inversion model of synthetic RMT TE-mode data. (d) Single inversion result of synthetic ERT data. (e) Joint inversion model of synthetic RMT TE-mode and ERT data. (f) Table showing RMS. The joint inversion model fits both data sets reasonably well and shows the highest level of detail among all the inversions. A synthetic model (Fig. 11b) was designed based on these interpretations to evaluate the suggested resolution properties of the joint inversion model. Synthetic lake-floor ERT and RMT data sets were generated at the same observation positions as the field data and contaminated with noise corresponding to the data error floors used in the inversions of the field data. Furthermore, the same data editing was used in subsequent inversions to better understand the inversions of field data. The results of single inversions (Figs 11c and d) and joint inversion (Fig. 11e) of the synthetic data reflect well how the field data were modeled in the single and joint inversions. The RMT single inversion resolved the facture zone at 600–700 m distance. Because of the low data coverage and structural complexity of the fracture zones at 500–700 m distance, the ERT single inversion has only incorrectly predicted one fracture zone (Fig. 11d). The joint inversion, however, resolves both fracture zones (Fig. 11e). The approximate maximal penetration depth of the combined data sets is also presented in Fig. 11(e) to evaluate the joint inversion results. 5 DISCUSSION AND CONCLUSIONS 5.1 Lake-floor ERT and boat-towed RMT joint inversion and useful techniques For the first time, a joint inversion of boat-towed RMT and lake-floor ERT data has been implemented based on a few modifications of an existing joint inversion code and a well-designed field case study. The implementation was tested with both synthetic and field data sets. All the results demonstrated that the joint inversion of boat-towed RMT TE-mode and lake-floor ERT data performs better than individual inversions. This is mainly due to: (1) lake-floor ERT and boat-towed RMT data together lead to improved data coverage of geological targets. (2) The sensitivities of the two methods complement each other (Candansayar & Tezkan 2008; Kalscheuer et al.2010). (3) The two methods are affected differently by noise. This also explains why a joint inversion model can fit both data sets, but the model inverted from one data set cannot explain the other data set. The inversion code EMILIA offers two different weighting schemes. One allows modifying the weights of different data sets and the other one allows modifying the weights on smoothness constraints in horizontal and vertical directions. Stronger weight on the data set, which has the lower number of data points, combines information from both data sets more evenly and results in a model that can explain both data sets (Fig. 7b). Weighting of smoothness constraints in different model directions is often required to explain data properly and, thus, it may eliminate artefacts produced by inversion with equal weights. In this study, this technique was used to remove possible artefacts in the vertical direction of the model (Fig. 7c). However, this did not influence the conductive zones at 550–700 m distance along the profile, i.e. the fracture system NE-1 and the possible fracture system EW-5 (Fig. 7c), which were already observed in the model in Fig. 7(b). This supports the hypothesis that the fracture zones NE-1 and EW-5 are resolved. A priori information provides known model structure to the inversion, such that ambiguity can be reduced in those parts of the model that are related to the a priori information (Siripunvaraporn et al.2005). This may lead to more accurate estimates of the model parameters that were determined with higher uncertainty without using a priori information. In the joint inversion using a priori information, the water resistivities, obtained from direct measurements, were fixed. Additionally, the smoothness constraints were decoupled along the lower boundary of the water body in order to correctly estimate the resistivity and thickness of the conductive lake sediments. All those setups together improved the inversion results for the fracture systems (Fig. 7d). 5.2 Evaluation of joint inversion and interpretation Exploration depth investigations, data misfit analyses and model resolution and error analyses were our main methods to evaluate the inversion models which had a reasonable RMS. Using only RMS as a standard evaluation tool is far from being sufficient. The well-constrained part of an inversion model is distinguished by being located within the depth of penetration range and represented by data with low misfit, focused model resolution and reasonably low model errors (Figs 7d, 8f and 9f and Tables 1 and 2). In this way, the well-constrained part of the model is deemed to be representative of the subsurface structures. One fracture zone (NE-1) is well delineated in this study (fracture zones EW-5 and EW-7 are also inferred although with a lower level of confidence) using the joint inversion approach applied to the boat-towed RMT and lake-floor ERT data sets. A synthetic test based on the interpretation results and another test just incorporating land RMT data (Fig. 12) proved that the joint inversion of the boat-towed RMT and lake-floor ERT data was indeed capable of resolving fractures in the bedrock. Detailed information about the geometry of three fracture zones (NE-1, EW-5 and EW-7) together with thick sediments at 200–300 m along the profile are new useful findings for SKB to consider in their ongoing research and can be drilled to be confirmed. However, the fracture zones NE-3 and NE-4 below the deep water and thick sediments cannot be resolved, because resolution underneath the conductive water and more conductive sediments is dramatically reduced. The promising results from the case study indicate that using boat-towed RMT and lake-floor ERT together would have higher resolution, if water and sediments were less conductive. Figure 12. View largeDownload slide Comparison between the joint inversion models inverted with boat-towed RMT data and without boat-towed RMT data. (a) Joint inversion with boat-towed RMT data, the figure is also shown in Fig. 6(c). (b) Joint inversion without boat-towed RMT data. (c) Table showing RMS. Both inversions have the same setup with regard to mesh, weighting and smoothness constraints. Apparently, boat-towed RMT data contribute information to the joint inversion model underneath the shallow water parts. Figure 12. View largeDownload slide Comparison between the joint inversion models inverted with boat-towed RMT data and without boat-towed RMT data. (a) Joint inversion with boat-towed RMT data, the figure is also shown in Fig. 6(c). (b) Joint inversion without boat-towed RMT data. (c) Table showing RMS. Both inversions have the same setup with regard to mesh, weighting and smoothness constraints. Apparently, boat-towed RMT data contribute information to the joint inversion model underneath the shallow water parts. Figure B1. View largeDownload slide (a) Inversion model computed from synthetic RMT determinant impedance data. (b) Inversion model computed from joint inversion of synthetic RMT determinant impedance and ERT data. (c) Inversion model computed from synthetic RMT TM-mode data. (d) Joint inversion model for synthetic RMT TM-mode and ERT data. (e) Inversion model for RMT TE- and TM-mode data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. The true model for which the synthetic data were computed is shown in Fig. 2(a). Symbols are the same as in Fig, 2. Figure B1. View largeDownload slide (a) Inversion model computed from synthetic RMT determinant impedance data. (b) Inversion model computed from joint inversion of synthetic RMT determinant impedance and ERT data. (c) Inversion model computed from synthetic RMT TM-mode data. (d) Joint inversion model for synthetic RMT TM-mode and ERT data. (e) Inversion model for RMT TE- and TM-mode data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. The true model for which the synthetic data were computed is shown in Fig. 2(a). Symbols are the same as in Fig, 2. 5.3 Prospect for lake-floor ERT and boat-towed RMT joint inversion With increasing population, higher demands will be imposed on future transportation systems not just on land but also under lakes and rivers. In Sweden, an area of about 40 000 km2 is covered by inland water, constituting around 95 700 lakes. Inevitably, parts of the transportation system have to pass either over or under the water. In the planning phase, site investigations of these areas will benefit from boat-towed RMT to design the optimal profile for lake-floor ERT, because boat-towed RMT is much more efficient than ERT. Then, the boat-towed RMT and lake-floor ERT joint inversion can further be used to map the most interesting area. This joint approach is more capable of resolving weakness zones effectively, for examples, unconsolidated sediments, faults and fracture zones, than single-method investigations based on either ERT or RMT. Thus, joint inversion of lake-floor ERT and boat-towed RMT data is a better way to improve model resolution in investigations of geological structure below shallow water bodies with those two geophysical methods at the present stage. This method can effectively resolve underwater structures and reduce the cost of pre-investigations in underwater infrastructure planning. For deeper (>10 m depending on the case) and conductive (<5 Ωm) water bodies, the RMT method will not provide sufficient depth penetration to explore structures below the water bodies. In such cases, a combination of boat-towed controlled-source RMT (in the frequency range of 1–10 kHz, Wang et al.2017), or boat-towed transient EM with a transmitter loop of adequate size (e.g. Barrett et al.2005; Hatch et al.2010; Mollidor et al.2013; Bekesi et al.2014) and lake-floor ERT methods should be employed. Acknowledgements This work was carried out within the frame of TRUST2.2 (http://trust-geoinfra.se) project sponsored by Formas (project number: 25220121907), SGU, BeFo, SBUF/Skanska, FQM and NGI. The ERT data were acquired within the TRUST 2.1 and 4.2 (http://trust-geoinfra.se) projects funded by Formas, BeFo and SBUF/Skanska. SKB provided access to the site and provided support for the measurements. Nova FoU is thanked for partly funding this work. UPPMAX (project numbers: SNIC 2016-1-38 and SNIC 2016-1-439) is acknowledged for providing computation facilities. M. Ronczka and P. Olsson are thanked for providing the ERT data. A.A. Pfaffhuber, M. Everett and one anonymous reviewer are thanked for their constructive and useful comments. S. Wang thanks the China Scholarship Council (201306370039), Formas (25220121907), National Natural Science Foundation of China (41227803) and National High Technology Research and Development Program of China (2012AA09A20105) for supporting his PhD study at Uppsala University. REFERENCES Abubakar A., Gao G., Habashy T.M., Liu J., 2012. Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data, Inverse Prob. , 28( 5), 055016, doi:10.1088/0266-5611/28/5/055016. Google Scholar CrossRef Search ADS   Almén K.E., Stenberg L., 2005. Characterisation Methods and Instruments . Experiences from the construction phase (no. SKB-TR-05-11). Äspö Hard Rock Laboratory. Swedish Nuclear Fuel and Waste Management Company. Barrett B., Heinson G., Hatch M., Telfer A., 2005. River sediment salt-load detection using a water-borne transient electromagnetic system, J. Appl. Geophys. , 58( 1), 29– 44. Google Scholar CrossRef Search ADS   Bastani M., 2001. EnviroMT—a new controlled source/radio magnetotelluric system, PhD thesis , Uppsala University. Bastani M., Hübert J., Kalscheuer T., Pedersen L.B., Godio A., Bernard J., 2012. 2D joint inversion of RMT and ERT data versus individual 3D inversion of full tensor RMT data: an example form Trecate site in Italy, Geophysics , 77( 4), WB233– WB243. Google Scholar CrossRef Search ADS   Bastani M., Pedersen L.B., 2001. Estimation of magnetotelluric transfer functions from radio transmitters, Geophysics , 66( 4), 1038– 1051. Google Scholar CrossRef Search ADS   Bastani M., Persson L., Mehta S., Malehmir A., 2015. Boat-towed radio-magnetotellurics—a new technique and case study from the city of Stockholm, Geophysics , 80( 6), B193– B202. Google Scholar CrossRef Search ADS   Bekesi G., Telfer A., Woods J., Forward P., Burnell R., Hatch M., 2014. Quantitative measure of salt interception using in-river transient electromagnetic geophysics, Aust. J. Water Res. , 18( 1), 55– 69. Berglund J., Curtis P., Eliasson T., Olsson T., Starzec P., Tullborg E.L., 2003 Update of the Geological Model 2002 . Äspö Hard Rock Laboratory. Swedish Nuclear Fuel and Waste Management Company. Brodic B., Malehmir A., Juhlin C., Dynesius L., Bastani M., Palm H., 2015. Multicomponent broadband digital-based seismic landstreamer for near-surface applications, J. Appl. Geophys. , 123, 227– 241. Google Scholar CrossRef Search ADS   Brodic B., Malehmir A., Juhlin C., 2017. Delineating fracture zones using surface-tunnel-surface seismic data, P-S and S-P mode conversions, J. Geophys. Res. , doi:10.1002/2017JB014304. Candansayar M.E., Tezkan B., 2008. Two-dimensional joint inversion of radiomagnetotelluric and direct current resistivity data, Geophys. Prospect. , 56( 5), 737– 749. Google Scholar CrossRef Search ADS   Christensen N.B., 1998. The determination of electrical anisotropy using surface electric and electromagnetic methods, in 11th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems , doi:org/10.4133/1.2922495. Constable S.C., Parker R.L., Constable C.G., 1987. Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics , 52( 3), 289– 300. Google Scholar CrossRef Search ADS   Coutant O., Bernard M.L., Beauducel F., Nicollin F., Bouin M.P., Roussel S., 2012. Joint inversion of P-wave velocity and density, application to La Soufrière of Guadeloupe hydrothermal system, Geophys. J. Int. , 191( 2), 723– 742. Google Scholar CrossRef Search ADS   Cosma C., Olsson O., Keskinen J., Heikkinen P., 2001. Seismic characterization of fracturing at the Äspö Hard Rock Laboratory, Sweden, from the kilometer scale to the meter scale, Int. J. Rock Mech. Min. Sci. , 38( 6), 859– 865. Google Scholar CrossRef Search ADS   Dahlin T., Zhou B., 2006. Multiple-gradient array measurements for multichannel 2D resistivity imaging, Near Surf. Geophys ., 4( 2), 113– 123. Dahlin T., Loke M.H., Siikanen J., Höök M., 2014. Underwater ERT survey for site investigation for a new line for Stockholm Metro, in Near Surface Geoscience 2014–20th European Meeting of Environmental and Engineering Geophysics , doi:10.3997/2214-4609.20142060. Dey A., Morrison H.F., 1979. Resistivity modelling for arbitrarily shaped two-dimensional structures, Geophys. Prospect. , 27( 1), 106– 136. Google Scholar CrossRef Search ADS   Daily W., Ramirez A., Binley A., LaBrecque D., 2005. Electrical resistance tomography—theory and practice, in Near-Surface Geophysics , pp. 525– 550, ed. Butler D.K., Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Edwards L.S., 1977. A modified pseudosection for resistivity and IP, Geophysics , 42( 5), 1020– 1036. Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near‐surface materials by joint 2D inversion of dc resistivity and seismic data, Geophys. Res. Lett. , 30( 13), doi:10.1029/2003GL017370. Gallardo L.A., Meju M.A., 2004. Joint two‐dimensional DC resistivity and seismic travel time inversion with cross‐gradients constraints, J. Geophys. Res. , 109(B3), doi:10.1029/2003JB002716. Gallardo L.A., Meju M.A., 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification, Geophys. J. Int. , 169( 3), 1261– 1272. Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2011. Structure‐coupled multiphysics imaging in geophysical sciences, Rev. Geophys. , 49( 1), doi:10.1029/2010RG000330. Haber E., Gazit M.H., 2013. Model fusion and joint inversion, Surv. Geophys. , 34( 5), 675– 695. Google Scholar CrossRef Search ADS   Hatch M., Munday T., Heinson G., 2010. A comparative study of in-river geophysical techniques to define variations in riverbed salt load and aid managing river salinization, Geophysics , 75( 4), WA135– WA147. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints, Geophysics , 74( 6), R99– R109. Google Scholar CrossRef Search ADS   Huang H., 2005. Depth of investigation for small broadband electromagnetic sensors, Geophysics , 70( 6), G135– G142. Google Scholar CrossRef Search ADS   Jones A.G., Groom R.W., 1993. Strike-angle determination from the magnetotelluric impedance tensor in the presence of noise and local distortion: rotate at your peril!, Geophys. J. Int. , 113( 2), 524– 534. Google Scholar CrossRef Search ADS   Kalscheuer T., Pedersen L.B., 2007. A non-linear truncated SVD variance and resolution analysis of two-dimensional magnetotelluric models, Geophys. J. Int. , 169( 2), 435– 447. Google Scholar CrossRef Search ADS   Kalscheuer T., Pedersen L.B., Siripunvaraporn W., 2008. Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents, Geophys. J. Int. , 175( 2), 486– 514. Google Scholar CrossRef Search ADS   Kalscheuer T., García Juanatey M.A., Meqbel N., Pedersen L.B., 2010. Non-linear model error and resolution properties from two-dimensional single and joint inversions of direct current resistivity and radiomagnetotelluric data, Geophys. J. Int. , 182( 3), 1174– 1188. Google Scholar CrossRef Search ADS   Kalscheuer T., Bastani M., Donohue S., Persson L., Pfaffhuber A.A., Reiser F., Ren Z.Y., 2013. Delineation of a quick clay zone at Smørgrav, Norway, with electromagnetic methods under geotechnical constraints, J. Appl. Geophys. , 92, 121– 136. Google Scholar CrossRef Search ADS   Kalscheuer T. et al., 2015. Joint inversions of three types of electromagnetic data explicitly constrained by seismic observations: results from the central Okavango Delta, Botswana. Geophys. J. Int. , 202( 3), 1429– 1452. Google Scholar CrossRef Search ADS   Lines L.R., Schultz A.K., Treitel S., 1988. Cooperative inversion of geophysical data, Geophysics , 53( 1), 8– 20. Google Scholar CrossRef Search ADS   Loke M.H., 1999. Electrical imaging surveys for environmental and engineering studies, in A Practical Guide to 2-D and 3-D Surveys . Loke M.H., 2002. Rapid 2D resistivity forward modelling using the finite-difference and finite-element methods. Loke M.H., Lane J.W., 2004. Inversion of data from electrical resistivity imaging surveys in water-covered areas, Explor. Geophys. , 35, 266– 271. Google Scholar CrossRef Search ADS   Makurat A., Løset F., Hagen A.W., Tunbridge L., Kveldsvik V., Grimstad E., 2006. A Descriptive Rock Mechanics Model for the 380–500 m Level . Report R-02-11, SKB. Malehmir A. et al., 2015. Planning of urban underground infrastructure using a broadband seismic landstreamer—tomography results and uncertainty quantifications from a case study in southwestern Sweden, GeoPhysics , 80( 6), B177– B192. Google Scholar CrossRef Search ADS   Mehta S., Bastani M., Malehmir A., Pedersen L.B., 2017. Resolution and sensitivity of boat-towed RMT data to delineate fracture zones–example of the Stockholm bypass multi-lane tunnel, J. Appl. Geophys. , 139, 131– 143. Google Scholar CrossRef Search ADS   Meju M.A., 1996. Joint inversion of TEM and distorted MT soundings: some effective practical considerations, Geophysics , 61( 1), 56– 65. Google Scholar CrossRef Search ADS   Menke W., 1989. Geophysical Data Analysis: Discrete Inverse Theory , International Geophysics Series , Vol. 45. Academic Press, London. Mollidor L., Tezkan B., Bergers R., Löhken J., 2013. Float‐transient electromagnetic method: in‐loop transient electromagnetic measurements on Lake Holzmaar, Germany, Geophys. Prospect. , 61( 5), 1056– 1064. Google Scholar CrossRef Search ADS   Monteiro Santos F.A., Andrade Afonso A.R., Dupis A., 2006. 2D joint inversion of dc and scalar audio-magnetotelluric data in the evaluation of low enthalpy geothermal fields, J. Geophys. Eng. , 4( 1), 53, doi:10.1088/1742-2132/4/1/007. Google Scholar CrossRef Search ADS   Moorkamp M., Heincke B., Jegen M., Roberts A.W., Hobbs R.W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data, Geophys. J. Int. , 184( 1), 477– 493. Google Scholar CrossRef Search ADS   Moorkamp M., Jones A.G., Fishwick S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data, J. Geophys. Res. , 115( B4), B04318, doi:10.1029/2009JB006369. Google Scholar CrossRef Search ADS   Moorkamp M., Lelievre P.G., Linde N., Khan A., 2016. Integrated Imaging of the Earth . Wiley Press. Google Scholar CrossRef Search ADS   Pedersen L.B., Engels M., 2005. Routine 2D inversion of magnetotelluric data using the determinant of the impedance tensor, Geophysics , 70( 2), G33– G41. Google Scholar CrossRef Search ADS   Ronczka M., Hellman K., Günther T., Wisén R., Dahlin T., 2017. Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory, Solid Earth , 8, 671– 682. Google Scholar CrossRef Search ADS   Rhén I., Gustafsson G., Stanfors R., Wikberg P., 1997. Äspö HRL-Geoscientific Evaluation 1997/5. Models Based on Site Characterization 1986–1995  (no. SKB-TR–97-06). Swedish Nuclear Fuel and Waste Management Company. Sasaki Y., 1989. Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data, Geophysics , 54( 2), 254– 262. Google Scholar CrossRef Search ADS   Schlumberger M., 1939. The application of telluric currents to surface prospecting, EOS, Trans. Am. Geophys. Un. , 20( 3), 271– 277. Google Scholar CrossRef Search ADS   Silvester P.P., Ferrari R.L., 1996. Finite Elements for Electrical Engineers . Cambridge University Press. Google Scholar CrossRef Search ADS   Stanfors R., Rhén I., Tullborg E.L., Wikberg P., 1999. Overview of geological and hydrogeological conditions of the Äspö hard rock laboratory site, Appl. Geochem. , 14( 7), 819– 834. Google Scholar CrossRef Search ADS   Sudha A., Tezkan B., Siemon B., 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter‐borne electromagnetic data, Geophys. Prospect. , 62( 3), 597– 614. Google Scholar CrossRef Search ADS   Siripunvaraporn W., Egbert G., 2000. An efficient data-subspace inversion method for 2D magnetotelluric data, Geophysics , 65( 3), 791– 803. Google Scholar CrossRef Search ADS   Siripunvaraporn W., Egbert G., Lenbury Y., Uyeshima M., 2005. Three-dimensional magnetotelluric inversion: data-space method, Phys. Earth Planet. Inter. , 150( 1), 3– 14. Google Scholar CrossRef Search ADS   Spies B.R., 1989. Depth of investigation in electromagnetic sounding methods, Geophysics , 54( 7), 872– 888. Google Scholar CrossRef Search ADS   Swift C.M. Jr, 1967. A magnetotelluric investigation of an electrical conductivity anomaly in the South Western United States, PhD thesis , MIT, Cambridge press. Tezkan B., Goldman M., Greinwald S., Hördt A., Müller I., Neubauer F.M., Zacher G., 1996. A joint application of radiomagnetotellurics and transient electromagnetics to the investigation of a waste deposit in Cologne (Germany), J. Appl. Geophys. , 34( 3), 199– 212. Google Scholar CrossRef Search ADS   Tezkan B., Hördt A., Gobashy M., 2000. Two-dimensional radiomagnetotelluric investigation of industrial and domestic waste sites in Germany, J. Appl. Geophys. , 44( 2), 237– 256. Google Scholar CrossRef Search ADS   TRUST, 2013. Available at: http://trust-geoinfra.se/, (accessed 2016 July 20). Tryggvason A., Linde N., 2006. Local earthquake (LE) tomography with joint inversion for P‐and S‐wave velocities using structural constraints, Geophys. Res. Lett. , 33( 7), L07303, doi:10.1029/2005GL025485. Google Scholar CrossRef Search ADS   Turberg P., Müller I., Flury F., 1994. Hydrogeological investigation of porous environments by radio magnetotelluric-resistivity (RMT-R 12–240 kHz), J. Appl. Geophys. , 31( 1-4), 133– 143. Google Scholar CrossRef Search ADS   Vozoff K., Jupp D.L.B., 1975. Joint inversion of geophysical data, Geophys. J. R. Astr. Soc. , 42, 977– 991. Google Scholar CrossRef Search ADS   Wang S., Bastani M., Kalscheuer T., Malehmir A., Dynesius L., 2017. Controlled source boat-towed radio-magnetotellurics for site investigation at Äspö Hard Rock Laboratory, Southeastern Sweden, in 79th EAGE Conference and Exhibition 2017 , doi:10.3997/2214-4609.201700565. Wikberg P., Gustafson G., Rhén I., Stanfors R., 1991. Evaluation and Conceptual Modelling Based on the Pre-investigations 1986–1990  (no. SKB-TR–91-22). Aespoe Hard Rock Laboratory.  Swedish Nuclear Fuel and Waste Management Company. Xiong B., Wang S.G., 2011. Wave-number domain features of primary field of H− Hz arrangement wide field electromagnetic method, in International Conference on Instrumentation, Measurement, Circuits and Systems (ICIMCS 2011) . ASME Press, doi:10.1115/1.859902.paper196. Xu S.Z., Duan B.C., Zhang D.H., 2000. Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2.5 D electrical modelling, Geophys. Prospect. , 48( 5), 789– 796. Google Scholar CrossRef Search ADS   Yogeshwar P., Tezkan B., Israil M., Candansayar M.E., 2012. Groundwater contamination in the Roorkee area, India: 2D joint inversion of radiomagnetotelluric and direct current resistivity data, J. Appl. Geophys. , 76, 127– 135. Google Scholar CrossRef Search ADS   Zhang P., Roberts R.G., Pedersen L.B., 1987. Magnetotelluric strike rules, Geophysics , 52( 3), 267– 278. Google Scholar CrossRef Search ADS   Zonge K., Wynn J., Urquhart S., 2005. Resistivity, induced polarization, and complex resistivity, in Near-Surface Geophysics , pp. 265– 300, ed. Butler, D.K., Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   APPENDIX A: SENSITIVITY CALCULATION Two types of ERT field data are usually collected in field applications. One is apparent resistivity (eq. 7), which is an average resistivity of the part of the subsurface penetrated by the ERT current system and equal to the real material resistivity when the subsurface is a homogeneous half-space. Another one is resistance (eq. 8). Both of them can be used to invert for a resistivity model. In EMILIA, apparent resistivities are transformed to logarithmic values to stabilize and speed up the inversion. Similarly, model resistivities are transformed to logarithmic values to speed up the inversion and avoid negative resistivity values. Thus, sensitivities are computed for logarithmic apparent resistivities with regard to logarithmic model resistivities yielding the following relationship to the sensitivities of the original quantities:   $$\frac{\partial \ \log _{10} (\rho _{a})}{\partial \ \log _{10} (\rho)} = - \frac{\partial \ \log _{10} (\rho _{a})}{\partial \ \log _{10} (\sigma )} = - \frac{\sigma }{\rho _{a}} \frac{\partial \rho _{a}}{\partial \sigma },$$ (A1)where σ = 1/ρ is conductivity. However, when resistances are used in inversion, it is less meaningful than for apparent resistivities to use logarithmic values, because resistances could be negative. For resistance data, sensitivities are calculated with regard to logarithmic model resistivities establishing the following relationship to the sensitivities of the original quantities   \begin{eqnarray} \frac{{\partial R}}{{\partial \log_{10} (\rho)}} &=&- \frac{{\partial R}}{{\partial \log_{10} (\sigma)}} = - \frac{{\partial R}}{{\partial \left( {\frac{{\ln \sigma }}{{\ln 10}}} \right)}} = - \ln 10\frac{{\partial R}}{{\partial \ln \sigma }}\nonumber\\ &=& - \ln 10 \cdot \sigma \frac{{\partial R}}{{\partial \sigma }}. \end{eqnarray} (A2) EMILIA was modified to use both types of data in inversion. This function makes EMILIA more flexible to utilize field data. APPENDIX B: INVERSIONS OF TM-MODE DATA We have also investigated inversions of synthetic TM mode, TE and TM modes and determinant impedance data computed for the true model in Fig. 2(a). However, the inversion results are only shown in this appendix, because we focus on the results of the TE-mode inversions. In general, the TM mode and determinant impedance single inversions are not in as good structural agreement with the true model as the TE-mode single inversion model, especially the inversion result of the determinant impedances is inferior to that of the TE-mode impedances (Fig. B1a and c). For the joint inversions with the ERT data, good resolution can still be observed in Fig. B1b (joint inversion with RMT determinant impedance data), B1d (joint inversion with RMT TM-mode impedance data) and B1f (joint inversion with RMT TE- and TM-mode impedance data). However, these joint inversion results are not as good as the model from joint inversion of TE-mode and ERT data (Fig. 2d). This is mainly so, because the TM-mode and ERT data have similar directions of electrical current flow. Thus compared to a single inversion of ERT data, the inclusion of TM-mode data may not add significant new information. Since determinant impedance data are geometrical averages of TE- and TM-mode impedances, this effect is also observed in joint inversions of determinant impedance and ERT data to some extent. The joint inversions of the TE- and TM-mode data and of the TE-, TM-mode and ERT data show as good results as the joint inversion of the TE-mode and ERT data. Thus, in our case, the TM-mode data hardly contribute information to our joint inversion of RMT and ERT data (Fig. B1c) that is not already contained in the ERT data. Field data of the TM mode, TE and TM mode, and determinant impedances were also inverted in the case study (TM-mode and tipper data are shown in Fig. B2). In general, the models below the water surface from single inversions of TM mode and determinant impedances (Figs B3a and c) are not as good as that of the TE-mode impedances (Fig. 6a). For the joint inversions, two conductive zones can still be observed in Figs B3(b), (d) and (f). However, they are not as clear as in the model from joint inversion of TE-mode and ERT data (Fig 6c). Figure B2. View largeDownload slide Unedited raw data: (a) TM-mode apparent resistivity data. (b) TM-mode phase data. (c) Real part of vertical magnetic transfer function Ty (Hz/Hy, in our coordinate system). (d) Imaginary part of vertical magnetic transfer function Ty. Figure B2. View largeDownload slide Unedited raw data: (a) TM-mode apparent resistivity data. (b) TM-mode phase data. (c) Real part of vertical magnetic transfer function Ty (Hz/Hy, in our coordinate system). (d) Imaginary part of vertical magnetic transfer function Ty. Figure B3. View largeDownload slide (a) Inversion model computed from RMT determinant impedance field data. (b) Inversion model computed from joint inversion of RMT determinant impedance and ERT field data. (c) Inversion model computed from RMT TM-mode field data. (d) Joint inversion model for RMT TM-mode and ERT field data. (e) Inversion model for RMT TE and TM-mode field data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. For the joint inversion, RMS values of individual data sets are slightly higher than those of the single inversions. Symbols are the same as in Fig. 6. Figure B3. View largeDownload slide (a) Inversion model computed from RMT determinant impedance field data. (b) Inversion model computed from joint inversion of RMT determinant impedance and ERT field data. (c) Inversion model computed from RMT TM-mode field data. (d) Joint inversion model for RMT TM-mode and ERT field data. (e) Inversion model for RMT TE and TM-mode field data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. For the joint inversion, RMS values of individual data sets are slightly higher than those of the single inversions. Symbols are the same as in Fig. 6. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Joint inversion of lake-floor electrical resistivity tomography and boat-towed radio-magnetotelluric data illustrated on synthetic data and an application from the Äspö Hard Rock Laboratory site, Sweden

23 pages

/lp/ou_press/joint-inversion-of-lake-floor-electrical-resistivity-tomography-and-CLCOPYcG7D
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx414
Publisher site
See Article on Publisher Site

### Abstract

Summary The electrical resistivity tomography (ERT) method provides moderately good constraints for both conductive and resistive structures, while the radio-magnetotelluric (RMT) method is well suited to constrain conductive structures. Additionally, RMT and ERT data may have different target coverage and are differently affected by various types of noise. Hence, joint inversion of RMT and ERT data sets may provide a better constrained model as compared to individual inversions. In this study, joint inversion of boat-towed RMT and lake-floor ERT data has for the first time been formulated and implemented. The implementation was tested on both synthetic and field data sets incorporating RMT transverse electrical mode and ERT data. Results from synthetic data demonstrate that the joint inversion yields models with better resolution compared with individual inversions. A case study from an area adjacent to the Äspö Hard Rock Laboratory (HRL) in southeastern Sweden was used to demonstrate the implementation of the method. A 790-m-long profile comprising lake-floor ERT and boat-towed RMT data combined with partial land data was used for this purpose. Joint inversions with and without weighting (applied to different data sets, vertical and horizontal model smoothness) as well as constrained joint inversions incorporating bathymetry data and water resistivity measurements were performed. The resulting models delineate subsurface structures such as a major northeasterly directed fracture system, which is observed in the HRL facility underground and confirmed by boreholes. A previously uncertain weakness zone, likely a fracture system in the northern part of the profile, is inferred in this study. The fractures are highly saturated with saline water, which make them good targets of resistivity-based geophysical methods. Nevertheless, conductive sediments overlain by the lake water add further difficulty to resolve these deep fracture zones. Therefore, the joint inversion of RMT and ERT data particularly helps to improve the resolution of the resistivity models in areas where the profile traverses shallow water and land sections. Our modification of the joint inversion of RMT and ERT data improves the study of geological units underneath shallow water bodies where underground infrastructures are planned. Thus, it allows better planning and mitigating the risks and costs associated with conductive weakness zones. Electrical resistivity tomography (ERT), Radio-magnetotellurics (RMT), Joint inversion, Fractures, faults, and high strain deformation zones 1 INTRODUCTION Due to the non-linearity and ill-posedness of geophysical inverse problems, seeking a unique model based on geophysical data is a difficult task. Joint inversion is a useful technique to narrow down the range of possible models and has become increasingly popular for interpreting geophysical data sets related to both common and disparate material properties. Two strategies, direct parameter coupling and structural coupling, are mainly used for joint inversion with different types of geophysical model parameters (e.g. Moorkamp et al.2011; Abubakar et al.2012; Haber & Gazit 2013; Moorkamp et al.2016). Direct parameter coupled joint inversion is based on petrophysical parameter measurements in the laboratory and only valid in limited parts of the model domain (Lines et al.1988; Coutant et al.2012). Structurally coupled joint inversion using cross-gradient constraints (Gallardo & Meju, 2003, 2004, 2007, 2011; Tryggvason & Linde 2006; Hu et al.2009; Moorkamp et al.2010; Moorkamp et al.2016) is easy to implement and offers sufficient flexibility to reduce structural coupling where this is in contradiction to the data. Single-property joint inversion is naturally coupled by inverting data from different methods for the same type of material property, for example among different types of electromagnetic (EM) methods (Meju 1996; Kalscheuer et al.2015). For the joint inversion of radio-magnetotelluric (RMT) and electrical resistivity tomography (ERT) data, many relevant publications exist on algorithm developments and case studies. First, Vozoff & Jupp (1975) described the advantages of joint inversion of MT and ERT data by synthetic and field examples using 1-D models. Sasaki (1989) implemented a 2-D joint inversion of MT and ERT dipole–dipole data and also tested it on synthetic and field data sets. Monteiro Santos et al. (2006) evaluated a geothermal field using 2-D joint inversion of MT and ERT data by the technique proposed by Sasaki (1989). Candansayar & Tezkan (2008) developed an algorithm for 2-D joint inversion of RMT and ERT data, which improves the resolution in modeling a fracture zone both in synthetic and field cases. This algorithm was also used to evaluate the impact of sewage irrigation and groundwater contamination in India by Yogeshwar et al. (2012). Kalscheuer et al. (2010) analysed the resolution and variance properties of 2-D resistivity models derived from single and joint inversions of ERT and RMT measurements. Bastani et al. (2012) investigated an oil contamination field in Italy by joint 2-D inversion of ERT and RMT data and compared the results with 3-D single inversion of RMT data. By studying a quick-clay site in Norway, Kalscheuer et al. (2013) demonstrated that the detectability of quick-clay zones could be improved by jointly inverting ERT and RMT data. For a few reasons, joint inversion of RMT (also MT data in general) and ERT data sets can lead to inverse models that are better constrained than those from the inversion of individual data sets. Due to spatial limitations, the ERT measurements may only constrain the shallow part of the subsurface (upper 20–30 m depending on the case) using electrode arrays that are a couple of hundred metres long. In cases with longer electrode arrays (>1 km) and large electrode spacing, ERT penetration depth may be substantially larger, but resolution for the shallow subsurface may be reduced. Hence, RMT measurements with a fixed range of signal frequencies of 14–250 kHz may constrain either correspondingly deeper or shallower parts of the subsurface. Thus, combining ERT and RMT data sets can lead to improved depth coverage of geological targets (Candansayar & Tezkan 2008). Furthermore, the sensitivities of the different methods complement each other; for example, ERT can resolve both resistive and conductive structures moderately well, while RMT is superior at resolving conductive structures but has negligible sensitivity to the resistivity of thin resistive units (Vozoff & Jupp 1975; Candansayar & Tezkan 2008; Kalscheuer et al.2010). Also, measurements by the different methods are differently affected by various types of field noise and couple differently to subsurface structures in the survey area. Finally, joint inversion of RMT and ERT data may to some extent help to determine electrical anisotropy parameters (Christensen 1998). Owing to the improvements in model constraints offered by joint inversions of ERT and RMT data and the increasing number of geo-engineering applications to investigate rock mass underneath shallow water bodies, we extended the joint inversion algorithm by Kalscheuer et al. (2010) to invert boat-towed RMT and lake-floor ERT data. The results of synthetic tests and a real-field application are presented in this work. This is a unique implementation and application, although lake-floor ERT (Dahlin et al.2014; Loke & Lane 2004) and boat-towed RMT (Bastani et al.2015; Mehta et al.2017) experiments have separately been shown, and joint inversion of ERT and RMT data collected on land has been done before (Candansayar & Tezkan 2008; Yogeshwar et al.2012; Bastani et al.2012; Kalscheuer et al.2013). To further improve the estimated joint inversion model, weighting techniques for the different data sets and incorporation of bathymetric data and water resistivity measurements as a priori constraints were employed. We applied the method to a real data set collected in an experiment at the Äspö Hard Rock Laboratory (HRL) 30 km north of the city of Oskarshamn in southeastern Sweden. The target was mainly a dipping fracture system (consisting of several subsets) that was intersected by the access tunnel of the HRL facility and was speculated to reach the lake floor: no clear evidence for this hypothesis was available from previous geoscientific studies. 2 THEORY 2.1 Radio-magnetotelluric method RMT is one type of passive-source EM methods where the signal sources are distant radio transmitters operating in the frequency range of 14 kHz –1 MHz. At sufficiently long distances, the EM signals are considered as plane waves simplifying the estimation of the electrical resistivity of near-surface structures (Bastani 2001; Bastani et al.2012; Tezkan et al.1996, 2000; Turberg et al.1994). RMT data comprise measurements of three components of the magnetic field (Hx, Hy, Hz) and two horizontal components of the electric field (Ex, Ey). In the frequency domain, the electric and magnetic field components are related through the complex-valued and frequency-dependent impedance tensor Z given as   $$\left[ {\begin{array}{@{}*{1}{c}@{}} {{E_x}}\\ {{E_y}} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} {{Z_{xx}}}&{{Z_{xy}}}\\ {{Z_{yx}}}&{{Z_{yy}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{H_x}}\\ {{H_y}} \end{array}} \right],$$ (1)and the vertical magnetic field and horizontal magnetic fields are related through the complex-valued and frequency-dependent vertical magnetic transfer function T given as   $${H_z} = \left[ {\begin{array}{*{20}{c}} {{T_x}}&{{T_y}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{H_x}}\\ {{H_y}} \end{array}} \right],$$ (2)where x and y represent measurement directions in the Cartesian coordinate system. The skin depth   $${\delta = 503\sqrt {\rho /f} },$$ (3)where f is frequency (unit is Hz) and ρ is resistivity (unit is Ωm), is typically used to evaluate the maximal exploration depth at the given signal frequency as 1.5 δ (Spies, 1989). In a layered medium, the resistivity ρ is replaced by an effective resistivity $$\tilde{\rho }$$ (Spies 1989; Huang 2005). In 2-D media, the diagonal impedance tensor elements are zero (Zxx = Zyy = 0) given that x is the strike direction of the 2-D geological structure and y is the profile direction. The apparent resistivities (ρxy/yx) and phases (φxy/yx) defined by Zxy (transverse electrical, TE mode, with currents flowing in x-direction) and Zyx (transverse magnetic, TM mode, with currents flowing in y- and z-directions) are used to estimate the general resistivity distribution:   $${{\rho _{xy/yx}} = \frac{1}{{\mu_0 \omega }}}{\left| {{Z_{xy/yx}}} \right|^2},$$ (4)  $${\varphi _{xy/yx}} = {\tan ^{ - 1}}\left( {\frac{{{\mathop{\rm Im}\nolimits} \left( {{Z_{xy/yx}}} \right)}}{{{\mathop{\rm Re}\nolimits} \left( {{Z_{xy/yx}}} \right)}}} \right)$$ (5)where μ0 and ω are permeability of free air and angular frequency. In areas where the strike direction is not well defined, inversion of determinant impedance data is more suitable to suppress 3-D effects in 2-D models (Perdersen & Engels 2005) than inversion of biased TE- or TM-mode data. Here, the determinant impedance is defined as   $${Z_{{\rm{DET}}}} = \sqrt {{Z_{xx}}{Z_{yy}} - {Z_{xy}}{Z_{yx}}} .$$ (6) 2.2 Electrical resistivity tomography Use of ERT dates back as far as the beginning of geophysics and originally was used to distinguish oil-saturated from water-saturated formations in hydrocarbon exploration (Schlumberger 1939). In near-surface geophysics, the ERT method has been used for hydrogeological, mining and geohazard investigations (Daily et al.2005). In ERT measurements, a direct current (I) is injected into the ground by means of two electrodes (one can be far away from the measurement area), and a potential difference or voltage (U) owing to this current flow is measured between two potential electrodes. The Wenner, dipole–dipole, Schlumberger and gradient (Dahlin & Zhou 2006) electrode configurations are commonly used for ERT surveys (Zonge et al.2005). ERT data are presented as pseudo-sections of apparent resistivity or resistance, which are defined as   $${\rho _a} = k\frac{U}{I}$$ (7)or   $$R = \frac{U}{I},$$ (8)respectively. Here, k is a geometric factor, which depends on the array type and electrode spacing. These measurements are geometrically weighted averages of the true resistivity distribution of the subsurface. For current injection in a homogeneous Earth, the section of the earth above the median depth of investigation (Loke 1999) has the same influence on the measured potential as the section below. This determines roughly to which depth anomalous geological structures can be investigated using a given array type (Loke 1999), and the maximal depth of investigation differs from array to array. 2.3 Inverse theory We used the code Electro-Magnetic Inversion using Least Intricate Algorithms (EMILIA, Kalscheuer et al.2008, 2010) for our implementation, carrying out synthetic tests and inverting the field data. EMILIA uses the finite-difference method to solve the Poisson equation (Dey & Morrison 1979; Kalscheuer et al.2010) and the Helmholtz equations (Siripunvaraporn & Egbert 2000; Kalscheuer et al.2008) for the ERT and RMT methods, respectively. In an inverse problem, an objective function of the form (Menke 1989; Kalscheuer et al.2010):   \begin{eqnarray} \Phi \left( {{{\boldsymbol m}},{\lambda }} \right) &=& {\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{\boldsymbol m}} \right]} \right)^{T}}{{\boldsymbol W}}_{d}^{T}{{{\boldsymbol W}}_{d}}\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{\boldsymbol m}} \right]} \right) \nonumber\\ &&-\, Q_{d}^{\rm{*}}\,{ +\, \lambda }{\left( {{{\boldsymbol m}} - {{{\boldsymbol m}}^{r}}} \right)^{T}}{{\boldsymbol W}}_{m}^{T}{{{\boldsymbol W}}_{m}}\left( {{{\boldsymbol m}} - {{{\boldsymbol m}}^{r}}} \right),\end{eqnarray} (9)is to be minimized with regard to the set of model parameters contained in the model vector $${{\boldsymbol m}} = {\rm{\ }}{( {{m_1}, \ldots ,{\rm{\ }}{m_{m}}} )^{T}}$$. $${{\boldsymbol d}} = {\rm{\ }}{( {{d_1}, \ldots ,{\rm{\ }}{d_{N}}} )^{T}}$$ is a vector containing N field data points (including RMT and ERT data), and the vector $${{\boldsymbol F}}[ {{\boldsymbol m}} ]$$ contains the corresponding forward responses computed for a given model m. For RMT, the responses are apparent resistivity and phase, and for ERT, the response is apparent resistivity or resistance (for our application, we chose resistance as will be presented later). The superscript T denotes matrix transposition. $${{{\boldsymbol W}}_{d}} = {\rm{\ diag}}{( {\sigma _1^{ - 1}, \ldots ,{\rm{\ \ }}\sigma _{N}^{ - 1}} )^{T}}$$ is a data weighting matrix, where σi is the standard deviation of the observed data. $$\mathbf{W}_m^T\mathbf{W}_m = \alpha_y\mathbf{\boldsymbol\partial}_y^T\mathbf{\boldsymbol\partial}_y + \alpha_z\mathbf{\boldsymbol\partial}_z^T\mathbf{\boldsymbol\partial}_z$$ is the model regularization operator, which controls the simplicity of the inversion model $${{\boldsymbol m}}$$ and contains vertical and horizontal smoothness operators ∂y and ∂z, respectively, with manually adjustable weights αy and αz (Kalscheuer et al.2010). $${{{\boldsymbol m}}^{r}}$$ is a reference model, which is constructed from a priori information. The Lagrange multiplier λ balances data fit and model simplicity. $$Q_{d}^{\rm{*}}$$ is the target data misfit. We use the following definitions of data misfit Qd and root mean square error (RMS) (Kalscheuer et al.2013):   \begin{eqnarray}{Q_{d,sw}}\left[ {{\boldsymbol m}} \right] \!=\! \frac{{{N_d}}}{{\sum\nolimits_{j \!=\! 1}^{{N_{ds}}} {\sum\nolimits_{i \!=\! 1}^{{N_j}} {{{\left( {\frac{1}{{{w_{ji}}}}} \right)}^2}} } }}{\sum\nolimits_{j \!=\! 1}^{{N_{ds}}} {\sum\nolimits_{i \!=\! 1}^{{N_j}} {\left( {\frac{1}{{{w_{ji}}}}\frac{{{d_{ji}} \!-\! {F_{ji}}\left[ {{\boldsymbol m}} \right]}}{{{\sigma _{ji}}}}} \right)} } ^2}\nonumber\\ \end{eqnarray} (10)and   $${\rm{RMS}} = \sqrt {\frac{{{Q_{d,sw}}\left[ {{\boldsymbol m}} \right]}}{{{N_d}}}},$$ (11)where Nd is the total number of data points, Nds is the number of data sets, Nj is the total number of data points in data set j, wji is data set weighting factor and σji is the standard deviation of dji. By linearizing the forward operator in the vicinity of the model of iteration k, $${{{\boldsymbol m}}_{k}}$$, the original problem of minimizing the function $$\Phi ( {{{\boldsymbol m}},{\rm{\lambda }}} )$$ is simplified to minimizing a function $${\Phi ^{{\rm{quad}}}}( {{{\boldsymbol m}},{\rm{\lambda }}} )$$ which is quadratic in $${{{\boldsymbol m}}_{{k} + 1}}$$ (Menke 1989),   \begin{eqnarray} &&{ {\Phi ^{{\rm{quad}}}}\left( {{{{\boldsymbol m}}_{{k} + 1}},{\rm{\lambda }}} \right)}\nonumber\\ &&{\quad = \left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{{{\boldsymbol m}}_{k}}} \right] - {{\bf J}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}_{k}}} \right)} \right)^{T}}\nonumber\\ &&{\qquad \times \,{{\boldsymbol W}}_{d}^{T}{{{\boldsymbol W}}_{d}}\left( {{{\boldsymbol d}} - {{\boldsymbol F}}\left[ {{{{\boldsymbol m}}_{k}}} \right] - {{\bf J}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}_{k}}} \right)} \right)}\nonumber\\ &&{\qquad - Q_{d}^{\rm{*}}{\rm{\, +\, \lambda }}{\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}^{r}}} \right)^{T}}{{\boldsymbol W}}_{m}^{T}{{{\boldsymbol W}}_{m}}\left( {{{{\boldsymbol m}}_{{k} + 1}} - {{{\boldsymbol m}}^{r}}} \right) ,}\end{eqnarray} (12)where $${\boldsymbol{J}} = {\{ {\frac{{\partial {F_{\boldsymbol{i}}}[ {{{{\boldsymbol m}}_k}} ]}}{{\partial {m_j}}}} \}_{m = {m_k}}}$$ is the Jacobian matrix of partial derivatives consisting of N rows and M columns. Thus, i = 1, …, N and j = 1, …, M. In this work, the ERT computation of forward responses and sensitivities in EMILIA has been modified to allow for fields generated by subsurface electrodes. Two types of wavenumber selections for Fourier transformation along the strike direction and inverse Fourier transformation can be chosen in EMILIA (Xu et al.2000; Xiong & Wang 2011). $${{\boldsymbol F}}[ {{{{\boldsymbol m}}_{k}}} ]$$ has been changed to represent either resistances or apparent resistivities, because apparent resistivity is difficult to be defined meaningfully when electrodes are half on land and half under water. The corresponding formulations of sensitivities are presented in Appendix  A. The RMT modeling part in EMILIA did not require any modification for our study. Model error and resolution analysis is required to study model stability and closeness of an estimated model to the true model. The ideal model resolution matrix is an identity matrix of size M × M, which means all model parameters are perfectly resolved (Menke 1989). Generally, the estimated model parameters are weighted averages of the true ones. Thus, the row of the model resolution matrix that pertains to the cell under investigation has non-zero entries also off its diagonal entry (referred to as spread; Menke 1989). Typically, the entries of a row of the model resolution matrix are scaled by the respective cell dimensions yielding a so-called resolving kernel. 3 SYNTHETIC TEST 3.1 ERT forward response: comparison of EMILIA and RES2DMOD In order to evaluate our modification of the ERT part in EMILIA with electrodes located at arbitrary positions, we compared the forward responses computed by EMILIA and RES2DMOD (Loke 2002) for the case of underwater electrodes. The program RES2DMOD (Loke 1999) is a combination of the classic algorithms described by Dey & Morrison (1979) and Silvester & Ferrari (1996). The model used for comparison is originally from Loke (2002), which was slightly modified to achieve a better comparison. It is a two-layer model with the first layer having a resistivity of 50 Ωm (e.g. water) and a thickness of 100 m. The second layer (i.e. the confining half-space) has a resistivity of 100 Ωm. A 10 m × 10 m wide block with a resistivity of 2000 Ωm is located with its top at 5 m depth below the top of the confining half-space. In Fig. 1(a), only the second layer is shown for a good view of the anomaly in the model. Electrodes (red triangles) with a Wenner Alpha array configuration and a spacing of 1 m were placed at the interface, set as 0 m depth, between the first and second layers. The total length of the profile is 50 m. The Wenner Alpha array, which has the strongest signal strength among the common arrays, is normally used for field surveys and usually just referred to as the ‘Wenner’ array (Edwards 1977). Forward responses were computed using the same mesh in the central part of the model with the anomalous block and electrodes with both RES2DMOD and EMILIA (the discretization in the outer parts of the grid towards the boundary nodes is not shown in RES2DMOD). All the simulated data points (408 in total) show relative differences ((ρres2dmod − ρemilia)/ρres2dmod × 100 per cent) between the two codes of less than 1.8 per cent (Fig. 1b). Following the plotting convention for pseudo-sections of the Wenner Alpha array, the horizontal position and depth of each data point are calculated as the midpoint between the potential electrode and 0.519 times the electrode spacing, respectively (Edwards 1977). Figure 1. View largeDownload slide Synthetic model consisting of two conductive layers (interface at z = 0 m, only deeper layer shown) and a 2-D anomaly in the second layer (Fig. 1a) and relative response differences (Fig. 1b) between ERT Wenner forward responses generated by RES2DMOD and EMILIA for the model in (a). Red triangles indicate electrodes. In (b), the horizontal axis represents the middle position of the two potential electrodes and the vertical axis represents 0.519 times a, which is the distance between two abutting electrodes (Edwards 1977). Note that only over limited areas, such as the anomaly and interface between layers 1 and 2 at z = 0 m, relative differences are as high as 1.8 per cent. Figure 1. View largeDownload slide Synthetic model consisting of two conductive layers (interface at z = 0 m, only deeper layer shown) and a 2-D anomaly in the second layer (Fig. 1a) and relative response differences (Fig. 1b) between ERT Wenner forward responses generated by RES2DMOD and EMILIA for the model in (a). Red triangles indicate electrodes. In (b), the horizontal axis represents the middle position of the two potential electrodes and the vertical axis represents 0.519 times a, which is the distance between two abutting electrodes (Edwards 1977). Note that only over limited areas, such as the anomaly and interface between layers 1 and 2 at z = 0 m, relative differences are as high as 1.8 per cent. 3.2 Inversion: synthetic example for single and joint inversions Before presenting the use of the modified implementation for the field case study, a synthetic example is shown. A test model (Fig. 2a) was abstracted from the single inversion result of the ERT field data (shown later). In the model, very conductive sea water with a resistivity of 1.38 Ωm exists at shallow depth in the central part (160–600 m) of the profile overlying bedrock with a resistivity of 1000 Ωm. A fracture zone (100 Ωm) that extends from 600 to 650 m along the profile and dips to northwest at an angle of about 50° from the horizontal is contained in the model (since the code only can use rectangular mesh, the fracture zone has step-like boundaries). The positions of RMT and ERT stations are marked by red triangles and black dots, respectively (Figs 2b and c). On the lake, the RMT station spacing is around 10 m and, on land, it varies from 20 to 40 m (yielding 52 stations in total). However, the frequency range is fixed from 14 to 250 kHz (nine frequencies equally distributed in log space). ERT data were generated for a gradient array type. The spacing between two adjacent electrodes is 5 m. The total length of the ERT cable is 790 m (as for the field example, data from 156 electrodes were used and data from one electrode were discarded owing to coupling problems or instrumental malfunctioning). EMILIA was used to model the RMT responses on the surface including water and land data and the ERT responses on land and at the bottom of the water body. Then, we contaminated the synthetic data with Gaussian noise of 5 per cent on ERT resistance, 10 per cent on RMT apparent resistivity and 2.29° on RMT phase data, which are on levels comparable to those observed in our field data. Note that the noise levels selected for the RMT apparent resistivities and phases correspond to a 5 per cent noise level on the impedance tensor elements. This corresponds to 5 per cent on ERT resistance data, if we assume that the errors in the RMT measurements are mainly in the electrical field. An experiment was carried out to compare the differences between individual and joint inversions with the synthetic data sets. Figure 2. View largeDownload slide (a) Test model deduced from single inversion of ERT field data (see below). (b) Single inversion model of RMT TE-mode data. Red triangles represent RMT stations at the surface and white dashed line shows the approximate exploration depth of RMT signals. (c) Single inversion model of ERT data. Black dots represent ERT electrode positions and white dashed line shows the approximate exploration depth of ERT. (d) Joint inversion model for RMT TE-mode and ERT data. White dashed line shows the approximate joint exploration depth of RMT and ERT signals. (e)Table showing RMS misfits. The joint inversion model can fit both data sets equally well. Figure 2. View largeDownload slide (a) Test model deduced from single inversion of ERT field data (see below). (b) Single inversion model of RMT TE-mode data. Red triangles represent RMT stations at the surface and white dashed line shows the approximate exploration depth of RMT signals. (c) Single inversion model of ERT data. Black dots represent ERT electrode positions and white dashed line shows the approximate exploration depth of ERT. (d) Joint inversion model for RMT TE-mode and ERT data. White dashed line shows the approximate joint exploration depth of RMT and ERT signals. (e)Table showing RMS misfits. The joint inversion model can fit both data sets equally well. Only the joint inversions of RMT TE-mode and ERT data are shown in this paper. For the RMT TM mode and the ERT method, the directions of current flow in a 2-D model are relatively similar with current flowing exclusively and predominantly, respectively, in the plane of the profile (i.e. perpendicular to the strike direction in both horizontal and vertical directions). Thus, including RMT TM-mode and ERT data in a joint inversion may lead to inclusion of relatively similar model constraints. Hence, joint inversion of RMT TE-mode and ERT data, where the current flows are perpendicular and predominantly parallel (owing to the use of point sources, there is an additional component of ERT current flow along the strike direction) to the plane of the profile, respectively, leads to more complementary model constraints than joint inversion of RMT TM-mode and ERT data. However, RMT TM-mode and determinant data results are shown in Appendix  B together with a discussion of the preferred use of the TE-mode data in the joint inversions. The single inversion model of the RMT TE-mode data is shown in Fig. 2(b). The shape of the water body is resolved. However, the fracture zone is not clearly resolved. The individual inversion of ERT data gives a similar result (Fig. 2c). White dashed lines represent the maximal exploration depths of the RMT and ERT data sets in the models estimated using Spies’ (1989) method and Dahlin and Zhou's (2006) method, respectively. The maximal RMT exploration depth was estimated using the lowest signal frequency of 14 kHz and the vertical resistivity section underneath each station. However, rather than using 1.5 skin depths as proposed by Spies (1989) to estimate the depth to which structure can be detected, we employed one skin depth to retrieve a conservative estimate for the depth to which the model is well constrained by the data. The joint inversion of RMT TE-mode and ERT data generated a model very similar to the true model (Fig. 2d). Both water and fracture zones are more accurately reconstructed than in any individual inversions. Especially, the dipping angle of the fracture zone is clearly evident. The contrast between the dipping conductor at 150 m distance and its surroundings is not particularly well pronounced. We would not interpret it as a fracture zone, and this feature may mainly be caused by the noise we added to our synthetic data. 4 CASE STUDY In the frame of the TRUST (TRansparent Underground STructure; TRUST 2013) project, boat-towed RMT and lake-floor ERT data were collected. The main aim of the project is to develop and implement methods and techniques to obtain more accurate models for planning, designing and constructing urban underground transportation systems, such as tunnels, subways and trams (Bastani et al.2015; Brodic et al.2015; Malehmir et al.2015; Brodic et al.2017; Mehta et al.2017). Äspö was chosen as one of several common sites for testing the effectiveness of the combined use of geophysical, geochemical and hydrological methods. Boat-towed RMT (Bastani et al.2015) and lake-floor ERT surveys (Dahlin et al.2014), among other methods, were carried out in May 2015 with the objective of imaging geological structures below a lake individually and in combination. 4.1 Geological setting Äspö is located approximately 30 km north of the Oskarshamn archipelago near the shoreline of the Baltic Sea in southeastern Sweden (Fig. 3). Granitic rocks and diverse types of fractures or fracture zones dominate the site (Cosma et al.2001). Most fracture zones at Äspö are a result of reactivation of older structures and appear mainly brittle. The style of such fracture zone tends to depend on the nature of any older structure being reactivated (Stanfors et al.1999). The lake where we performed our measurements is connected to the Baltic Sea through a narrow water channel. In the northernmost part of the site, the Äspö HRL operated by SKB (Swedish Nuclear Fuel and Waste Management Company) is located. A nuclear power plant is about 500 m away from the surveyed line. An underground tunnel (partly shown as an NNW trending blue solid line in Fig. 3b) plunging at about 14 per cent over a length of about 1500 m connects the surface to the HRL as well as various smaller tunnels, ramps and one main shaft (Almén & Stenberg 2005). Figure 3. View largeDownload slide (a) Areal photo (Google map) and (b) geological map of the case study area (courtesy of SGU). RMT profiles are marked by different symbols with different colours. The RMT profile used for the joint inversion is marked by red circles on the lake and pink triangles on the land. ERT data were measured only for one profile, marked by green dots. Positions of fracture zones, which are documented by SKB, are shown in (b). However, they are projected to the ground surface based on limited tunnel and borehole observations combined with lost low-resolution geophysical data. The bedrock is mainly granitic. The general topography of the land part is approximately flat. Figure 3. View largeDownload slide (a) Areal photo (Google map) and (b) geological map of the case study area (courtesy of SGU). RMT profiles are marked by different symbols with different colours. The RMT profile used for the joint inversion is marked by red circles on the lake and pink triangles on the land. ERT data were measured only for one profile, marked by green dots. Positions of fracture zones, which are documented by SKB, are shown in (b). However, they are projected to the ground surface based on limited tunnel and borehole observations combined with lost low-resolution geophysical data. The bedrock is mainly granitic. The general topography of the land part is approximately flat. Based on the information obtained during tunnel construction, an NE-SW running fracture system (known as NE-1) exists below the northern end of the RMT lake profile (Fig. 3b). It can be regarded as a southern boundary to the laboratory and plays an important role for the geometries of the fractures and deformation zones that formed at the Äspö HRL site (Berglund et al.2003). At about 1300 m along the access tunnel, the NE-1 fracture system is intersected at approximately 180 m below ground surface (Almén & Stenberg 2005). Three main subsets comprise the fracture system and are about 60 m wide in total (Berglund et al.2003; Rhén et al.1997). The two southernmost subsets, trending NE and dipping to NW, can be described as highly fractured and hydraulically conductive zones (Stanfors et al.1999; Makurat et al.2006). The northern subset, approximately 28 m wide, is the most intensely fractured part of NE-1 and produced significant inflow of water to the tunnel during the construction phase (Stanfors et al.1999; Makurat et al.2006; Almén & Stenberg 2005). Diorite, fine-grained granite and greenstone are the main host rocks of the zone (Stanfors et al.1999; Berglund et al.2003). The central part of NE-1 is strongly altered to clay. Since NE-1 is strongly water bearing, various other fractures around the zone with various orientations (Fig. 3b) are also water bearing (Berglund et al.2003). The water is a mixture of non-saline and brackish sea water (Wikberg et al.1991). The dip of the NE-1 fracture zone in the tunnel intersection and adjusted by cored boreholes is about 65° with uncertainty, because information gaps exist between the tunnel and the boreholes (Berglund et al.2003). Further away from the tunnel and boreholes, the complex geological characteristics of the zone are only poorly determined. Furthermore, the NE-1 fracture zone is not exposed at the surface. Hence, geophysical information becomes especially important for resolving its extent and geometry (Berglund et al.2003). Three fracture zones namely EW-7, NE-4 and NE-3 in the southern half of the RMT profile below the lake (Fig. 3b) were indicated by refraction seismic (not publically available) and borehole data (Wikberg et al.1991; Stanfors et al.1999; Rhén et al.1997). The NE-4 fracture zone trends NE and consists of two continuous subsets with an approximate width of 40 m. In the access tunnel, the NE-3 fracture zone is 49 m wide. It dips steeply towards NNW (Stanfors et al.1999). Geological information suggests the existence of the fracture zone EW-5 (Fig. 3b, Wikberg et al.1991); however, it cannot be observed in the tunnel. In our study, new geophysical data along with bathymetry data and water resistivity measurements are used in 2D inversions to delineate the shape of the NE-1 fracture zone and others that are inferred such as EW-5. 4.2 Data acquisition and quality The instrument Enviro-MT (Bastani 2001) developed at Uppsala University (UU) was used for RMT data acquisition. This system was upgraded to a boat-towed data acquisition system in 2015 thanks to the collaboration between UU and the Geological Survey of Sweden (SGU, Bastani et al.2015). In our field measurements, the RMT station spacing was not fixed. On the lake the spacing was around 10 m, but on land the spacing varied from 20 to 40 m depending on where sufficient space was available to place the RMT system (Fig. 3a) resulting in 52 RMT stations in total. The signal frequency range of the received Very Low Frequency and Low Frequency radio transmitters is from 14 to 250 kHz, which as a result of further processing following Bastani & Pedersen (2001) gave nine frequencies equally distributed in log space. Compared to the land part, the RMT data acquisition on the lake was very efficient. Within two and a half hours, all the RMT stations on the lake were surveyed (Fig. 3a). However, two days were used to acquire the land RMT data (Fig. 3a). RMT stations surveyed and used for joint inversion are marked by red circles on the lake and pink triangles on land. The instrument Terrameter LS, designed by ABEM, was used to carry out ERT data acquisition using a multiple gradient array (Dahlin & Zhou 2006). The spacing between two adjacent electrodes was 5 m and the ERT electrode spread was 635 m long. The survey line was extended by roll-along to a total length of 790 m (after editing, data from 156 electrode positions out of a total of 157 electrode positions could be used). However, since the multielectrode cable was placed along the lake floor combined with some parts on land, the horizontal distance is 783.5 m along the profile. Along the lake floor, the electrode outlets of the cable were placed in the sediments (the electrode cable sank into the sediments when it was unrolled from a boat, because its density is higher than that of water). Since this provided sufficiently low coupling resistance, connecting electrodes to the cable was rendered unnecessary. The cable outlets on land were connected to plate electrodes, which were buried where glacial sediments were present and attached on the granite bedrock with a conductive gel on outcrops. More details of the ERT data acquisition can be found in Ronczka et al. (2017). ERT stations used for joint inversion are marked by green dots in Fig. 3(a). The field data sets are shown in Fig. 4. The RMT TE-mode apparent resistivity and phase data collected on water suggest that the sea water (characterized by the signals with frequencies >100 kHz) has even higher resistivity than the underlying sediments (characterized by the 50–100 kHz signals; Figs 4a and b), which together with the sea water strongly reduced the penetration depth of the RMT signals (eq. 3). The total thickness of water and sediments is small (<3.5 m) along 50 per cent of the part of the profile covered by water, for example, at 380–600 m distance along the profile (Fig. 4d), and around 20 per cent of the RMT stations are on land. Therefore, the RMT data can provide useful information for modeling the underlying resistive bedrock and the lake sediments. Our vertical magnetic transfer function data are not of sufficiently high quality to be used in the interpretation (convergence is bad when inverting magnetic transfer function data). Since the ERT cable was placed on the lake floor, the ERT penetration depth in the underlying resistive basement is not affected strongly by the overlying conductive water. The ERT data in Fig. 4(c) show low-resistance values in the centre of the line and high-resistance features at both ends of the line corresponding to outcropping bedrock. Before the data sets were inverted, noisy data were removed. Owing to the diffusive and integrative natures of EM and electrical fields, the criterion for removing data points was deviation from smooth response variation. For the inversion, a 10 per cent relative error floor was used for the RMT apparent resistivity data (30 per cent relative error floor was used for the land RMT apparent resistivity data to avoid influence from static shift), an absolute error floor of 2.29° was used for RMT phase data, and a 5 per cent relative error floor was used for the ERT data. Figure 4. View largeDownload slide Field measurements along the coincident RMT and ERT profiles: (a) RMT TE-mode apparent resistivity data, (b) RMT TE-mode phase data, (c) ERT resistance data and (d) elevation of lake floor. The lake-floor sediments have lower resistivity than water which reflects in RMT apparent resistivity (panel a) at frequencies of 50–100 kHz and at 300–500 m distances along the profile being lower than at the highest frequencies, which penetrate only into water. Figure 4. View largeDownload slide Field measurements along the coincident RMT and ERT profiles: (a) RMT TE-mode apparent resistivity data, (b) RMT TE-mode phase data, (c) ERT resistance data and (d) elevation of lake floor. The lake-floor sediments have lower resistivity than water which reflects in RMT apparent resistivity (panel a) at frequencies of 50–100 kHz and at 300–500 m distances along the profile being lower than at the highest frequencies, which penetrate only into water. 4.3 RMT strike analysis Based on the relevant geological information and previous studies of the site, the RMT data acquisition parameters (profile orientation and station spacing) were carefully chosen to resolve possible fracture zones (Fig. 3). However, dimensionality, distortion and strike analysis of the RMT data needs to be done to decompose the RMT field data into TE- and TM-mode data. A number of publications discussed the strike rules for the magnetotelluric impedance tensor in details (e.g. Zhang et al.1987; Jones & Groom 1993). In our study, the method proposed by Zhang et al. (1987) was used for dimensionality, distortion and strike analysis by calculating distortion parameters and strike directions from the impedance tensor and evaluating the fit of the distortion model to the field data. When applying this method to the RMT data, independent strike angles are obtained for each station and frequency. For our data, only one preferred direction of strike shows in a cumulative rose diagram (Fig. 5a), which is 83°E. Since we oriented the sensors of both the land and boat-towed RMT stations parallel and perpendicular to the profile direction, this strike direction is relative to the profile direction, and it corresponds to 75° East of geographic North. Swift (1967) skews of all of our RMT data (Fig. 5b, except for the impedance tensor of one station at the lowest signal frequency) are lower than 0.2 and most of them are lower than 0.1. Therefore, our RMT data can be considered to represent a 1-D/2-D structure and 3-D effects should not become problematic in 2-D inversions and in later interpretation of the 2-D resistivity model. According to Fig. 3(b), the RMT profile used for joint inversion is almost perpendicular to the speculated fracture zone orientation. Since the field measurements were well designed using the available fracture zone information and the joint inversion of RMT and ERT data requires coincident profile directions, it turned out that there was no need for rotation of the reference coordinate system to the strike direction. However, projection of the first five station locations onto a modeling profile perpendicular to the strike direction at the southeastern end was needed, because these stations are 35–50 m away from the profile (Fig. 3a). Figure 5. View largeDownload slide (a) Strike directions calculated with the method proposed by Zhang et al. (1987). (b) Swift skews. Both of them show the structure underneath the research area is 2-D with a preferred strike direction of 83° East with regard to the profile direction used in the field (corresponding to 75° East of geographic North). Figure 5. View largeDownload slide (a) Strike directions calculated with the method proposed by Zhang et al. (1987). (b) Swift skews. Both of them show the structure underneath the research area is 2-D with a preferred strike direction of 83° East with regard to the profile direction used in the field (corresponding to 75° East of geographic North). 4.4 Single and joint inversion results We carried out the inverse modeling in two steps. In the first step, we performed an Occam inversion using regular smoothness constraints (Constable et al.1987; Menke 1989) with a variable Lagrange multiplier and a 100 Ωm half-space initial model. In the second step, the final inversion model from step one was then used as a new initial model in an Occam inversion with additional Marquardt–Levenberg damping. Here, the Langrage multiplier for the smoothness constraint is kept fixed at the value that gave the lowest RMS in step one, while an optimal damping factor (values varied from 100.0 to 102.8 in our inversions) is determined in each iteration using a line search. Note, in a few cases using the average resistivity of the final inversion model from step one as a new initial model for step two resulted in a better data fit. Consequently, we considered this model for further considerations. The models from individual and joint inversions of the field data are shown in Fig. 6. Both individual inversions show a conductive water body, more conductive sediments saturated with saline water below the water body (ERT electrodes positions mark the water bottom in Fig. 6b) and a high-resistivity feature below the lake floor in the central part of the profile (at approximately 400 m along the profile) most possibly correlating with the small island east of the profile (Fig. 3). Besides, two low-resistivity zones (possibly fracture zones) appear at distances of 200–350 and 550–700 m along the profile in both single inversion models (Figs 6a and b). Comparing the two models (Figs 6a and b), the RMT data provide more information of the possible fracture zone at 550–600 m distance along the profile than the ERT data (white dashed lines indicate exploration depth in both models). However, there is a strong similarity between the conductive structures in both models. The joint inversion model (Fig. 6c) shows two relatively conductive zones at 550–700 m distance along the profile. Both of them are above the exploration depth (the white dashed line, which is the maximal exploration depth of the joint RMT and ERT data sets). The same consideration of exploration depth is suitable for the conductive zone at 200–350 m distance along the profile. The RMS misfits (in panel d of Fig. 6) of the joint and single inversion models are comparable; however, the joint inversion model can fit two different geophysical data sets simultaneously and more detailed structures are introduced in the model. This is because the null space (i.e. the part of the model space which is not determined by the data, e.g. Menke 1989) corresponding to the RMT TE data can partly be resolved by the ERT data and vice versa. This is also well proven by another experiment as follows. The models from single inversions of RMT and ERT data were used as initial models for single inversions of the ERT and RMT data, respectively, providing initial RMS misfits of 9 in either case. Considering the final RMS of 2.9 in the joint inversion (Fig. 6) implies that a single inversion model from one data set does not explain the other data set well, because each data set has its own limitations. In such a situation, joint inversion plays a key role to reduce the null space due to the differences in data coverage, sensitivity and noise influences. Figure 6. View largeDownload slide Inversion of field data: (a) Inversion model computed from RMT TE-mode data. (b) Inversion model computed from ERT data. (c) Joint inversion model for RMT TE-mode and ERT data. (d) Table showing RMS. Separate RMS values from joint inversion are slightly higher than those of the single inversions. However, the model fits both data sets with acceptable RMS. Figure 6. View largeDownload slide Inversion of field data: (a) Inversion model computed from RMT TE-mode data. (b) Inversion model computed from ERT data. (c) Joint inversion model for RMT TE-mode and ERT data. (d) Table showing RMS. Separate RMS values from joint inversion are slightly higher than those of the single inversions. However, the model fits both data sets with acceptable RMS. Although TE-mode data were chosen for this study, TM, TE and TM, and determinant data were also inverted for 2-D models. In Appendix  B, single inversion models of TM-mode and determinant impedances and joint inversions combined with ERT data are presented (see also Fig. B3). 4.5 Weighting of data sets and smoothness constraints ERT and RMT data sets usually have different data point densities and different sensitivities to model structures and are affected differently by 3-D anomalies implying the need of different weights on the ERT and RMT data sets (Kalscheuer et al.2013; Sudha et al. 2014). Since the ERT data set has more data points and better data coverage than the RMT data set, more weight on the ERT data set only results in the joint inversion model being close to the ERT single inversion model (Fig. 7a). The RMS values of the individual data sets also prove that the joint inversion is dominated by the ERT data (see the caption of Fig. 7). Moreover, higher weight on the RMT data set leads to inclusion of more information from the RMT data, which is superior to resolve conductors (in our case the possible fracture zone). Joint inversion with two times more weight on the RMT data set shows two conductive zones at 550–700 m distance along the profile (Fig. 7b). Those two conductive zones are consistent with the ones in the inversion (Fig. 6c) without any weighting. Other weights (1.5:1, 3:1 and 4:1) were also tried without giving better results than the one shown. Figure 7. View largeDownload slide (a) Joint inversion model using weight on ERT data two times higher than on RMT data. Apparently, this inversion is dominated by ERT data set. (b) Joint inversion model using weight on RMT data two times higher than on ERT data. Since the ERT data set has a higher number of data points, higher weight on the RMT data leads the joint inversion to adequately fit both data sets. (c) Joint inversion model when weight on horizontal model smoothness is twice the weight on vertical smoothness. (d) Joint inversion model constrained with bathymetry data and water resistivity measurement. (e) Zoom-in figure of part of the model in (d) marked by a rectangle. (f) Table showing RMS. All symbols in the subfigures are the same as in Fig. 6. All joint inversions are consistent with regard to the conductive zones along the profile except for the ERT dominated one. Figure 7. View largeDownload slide (a) Joint inversion model using weight on ERT data two times higher than on RMT data. Apparently, this inversion is dominated by ERT data set. (b) Joint inversion model using weight on RMT data two times higher than on ERT data. Since the ERT data set has a higher number of data points, higher weight on the RMT data leads the joint inversion to adequately fit both data sets. (c) Joint inversion model when weight on horizontal model smoothness is twice the weight on vertical smoothness. (d) Joint inversion model constrained with bathymetry data and water resistivity measurement. (e) Zoom-in figure of part of the model in (d) marked by a rectangle. (f) Table showing RMS. All symbols in the subfigures are the same as in Fig. 6. All joint inversions are consistent with regard to the conductive zones along the profile except for the ERT dominated one. Effects of different weights on the horizontal and vertical smoothing of the model (Kalscheuer et al.2010) were also evaluated using weights of one on both data sets. More weight on horizontal smoothing was used to generate structures that are extended predominantly in the horizontal direction. Twice more weight on horizontal smoothing than on the vertical smoothing did not smear out any of the vertical conductive zones at 550–700 m distance along the profile (Fig. 7c) suggesting that these conductive zones are required by the data. In this test with increased weight on horizontal smoothing, using the data set weights that were previously found to be optimal (i.e. two times more weight on the RMT data than on the ERT data) resulted in increased RMS misfits (total RMS = 3.61, RMT TE-mode data RMS = 3.05 and ERT data RMS = 4.32). This seems to be related to the fact that the RMT TE-mode data and ERT data have different sensitivities to horizontal and vertical resistivity contrasts thus necessitating adjustment of data set weights for changes in model smoothness constraints and vice versa. 4.6 Joint inversion using an a priori model During the RMT lake measurement, bathymetry data were also acquired. The water resistivity was measured at three different depths (Ronczka et al.2017). A three-layer model based on the water resistivity measurements (1.48 Ωm at 0.2 m depth, 1.37 Ωm at 3.8 m depth and 1.28 Ωm at 4.8 m depth.) can be used to represent the resistivity of the lake water. This information was used as additional constraints in the joint inversion by constructing an a priori model (mr in eq. 12). The a priori model had two structures, the water body with three layers and the rest of the model consisting of bedrock with a resistivity of 225 Ωm. The geometry of the low boundary of the water body in the a priori model was controlled by the bathymetry data. The bedrock resistivity in the a priori model was the average apparent resistivity value of presumed bedrock structure in the joint inversion models in Figs 6(c) and 7(a)–(c). Since the water resistivity was known at three different depths, it was fixed in the hope that this would constrain the remaining model parameters to find representative values. However, below the lake floor, more conductive sediments were observed in the inversion models without bathymetric constraints. Hence, in implementing the bathymetric constraints, the smoothness constraints along the lake floor were removed to guarantee the free change of the model parameters across the lake floor. By this way, the water resistivity and the lake floor depth were introduced in the inversion without influencing the parts of the model lacking the a priori information. The two-step inversion strategy presented in Section 4.4 was used. The initial model in step one is the a priori model built from the available information, and the a priori model and smoothness constraints constructed according to water resistivity, average apparent resistivity and bathymetry were used in both inversion steps. The joint inversion model of the RMT TE-mode and ERT data constrained by bathymetry and measurements of water resistivity (Fig. 7d) is similar to the results without these constraints. In this model, the shape of the lake-floor sediments is more focused (Fig. 7e). Also, the bedrock surface is better defined (as compared to the other joint inversions) and the resistivity values of bedrock are more representative of crystalline rocks (mostly above 10 000 Ωm). The conductive zones at 550–700 m distance along the profile are at least as well resolved as in the other joint inversions. In this inversion, the weights of the RMT TE-mode and ERT data sets and of the horizontal and vertical smoothness constraints were set to one. In alternative inversions using weights deviating from one (results not shown), the RMS misfits were higher. 4.7 Evaluation of inversion models We evaluate the joint inversion results by analysing (1) the exploration depth (already shown above), (2) the data misfit for each data point and (3) model resolution and error estimates. Data misfits ((dfield − demilia)/σfield, where σfield represents data uncertainties) are displayed for all the inversion models in Figs 8 and 9. In general, the RMT phase data have better fits than the RMT apparent resistivity data (Fig. 8). Compared with the RMT data misfit values (Fig. 8), the ERT data misfits are lower (Fig. 9). Importantly, joint inversion models without bathymetric and water resistivity constraints show data misfits comparable to those of the single inversion results (Figs 8 and 9). In particular for the RMT data, the joint inversion with bathymetric and water resistivity constraints shows higher RMS than the single inversions. However, a careful look at the data differences suggests that the higher frequency data corresponding to the water layer apparently show higher misfits (Fig. 8f). This is likely because the assumption of a three-layer water resistivity model is not accurate enough. A one-layer model was tested using an average of the water resistivities and a complex model was also tested using a linear decrease of the water resistivity with depth, but both inversions showed slightly degraded fits compared to the three-layer model. The conductive zones at 500–700 m distance are better resolved in the joint inversion model computed with constraints from bathymetry and water resistivity (Fig. 7d) yielding small data differences comparable to the single inversion results. This feature can be observed in Fig. 8(f) and Table 1, especially in the fit to the RMT apparent resistivity. Since the ERT data were collected on the lake sediments, the likely inaccurate assumption of water resistivity seems to have little influence and a careful look at the misfit reveals the improvement of data misfit (the elliptic zone marked in Fig. 9f as compared to Figs 9a–e). Thus, the joint inversion result constrained by bathymetry and water resistivity in Fig. 7(d) was chosen for interpretation. In general, the joint inversion improves the details of the inverted model and results in fitting both data sets with RMS values similar to those of single inversions. Figure 8. View largeDownload slide Normalised data differences between field and inverse modeled RMT data from (a) RMT TE data single inversion; (b) joint inversion without any weighting and a priori information; (c) joint inversion when weight on ERT data is two times higher than on RMT data; (d) joint inversion when weight on RMT data is two times higher than on ERT data; (e) joint inversion when weight in horizontal direction on the model is two times higher than in vertical direction and (f) joint inversion when constrained with bathymetry data and water resistivity measurements. Most joint inversions have RMS misfits comparable to the single inversion. In general, the phase is better fitted than the apparent resistivity. Crosses show the data points removed before the inversions were run (around 13 per cent of apparent resistivity data and 16 per cent of phase data). In the constrained joint inversion (f), a relatively large misfit is observed, because describing the water column with three layers is not accurate enough to represent the spatial variability of water resistivity. However, the part below the water zone gives comparable or even better [apparent resistivity in (f)] misfit compared to single inversions. The relevant RMS values are shown in Table 1. Figure 8. View largeDownload slide Normalised data differences between field and inverse modeled RMT data from (a) RMT TE data single inversion; (b) joint inversion without any weighting and a priori information; (c) joint inversion when weight on ERT data is two times higher than on RMT data; (d) joint inversion when weight on RMT data is two times higher than on ERT data; (e) joint inversion when weight in horizontal direction on the model is two times higher than in vertical direction and (f) joint inversion when constrained with bathymetry data and water resistivity measurements. Most joint inversions have RMS misfits comparable to the single inversion. In general, the phase is better fitted than the apparent resistivity. Crosses show the data points removed before the inversions were run (around 13 per cent of apparent resistivity data and 16 per cent of phase data). In the constrained joint inversion (f), a relatively large misfit is observed, because describing the water column with three layers is not accurate enough to represent the spatial variability of water resistivity. However, the part below the water zone gives comparable or even better [apparent resistivity in (f)] misfit compared to single inversions. The relevant RMS values are shown in Table 1. Figure 9. View largeDownload slide Normalised data differences between field and modeled ERT data for the same [except for the single inversion in (a)] inversion models as considered in Fig. 8. Joint inversions have RMS misfits comparable to the single inversion. The data marked by an ellipse at 600 m along the profile cause high misfit. Apparently, the joint inversion constrained by bathymetry and measurements of water resistivity shows the best fit over this area. Figure 9. View largeDownload slide Normalised data differences between field and modeled ERT data for the same [except for the single inversion in (a)] inversion models as considered in Fig. 8. Joint inversions have RMS misfits comparable to the single inversion. The data marked by an ellipse at 600 m along the profile cause high misfit. Apparently, the joint inversion constrained by bathymetry and measurements of water resistivity shows the best fit over this area. Table 1. Averages of data differences shown in Fig. 8. However, the three highest frequencies mainly relating to saline water are excluded in order to have a better comparison among the different types of inversion (cf. main text). The constrained joint inversion (Fig. 8f) has obviously the smallest data differences compared with the single inversion (Fig. 8a). Data differences  Fig. 8(a)  Fig. 8(b)  Fig. 8(c)  Fig. 8(d)  Fig. 8(e)  Fig. 8(f)  Average ρ difference  2.08  2.59  2.78  2.43  2.83  1.93  Average φ difference  0.79  1.16  1.91  0.94  1.26  1.24  Data differences  Fig. 8(a)  Fig. 8(b)  Fig. 8(c)  Fig. 8(d)  Fig. 8(e)  Fig. 8(f)  Average ρ difference  2.08  2.59  2.78  2.43  2.83  1.93  Average φ difference  0.79  1.16  1.91  0.94  1.26  1.24  View Large A model error and resolution analysis was then done for the joint inversion model constrained by bathymetry and water resistivity in Fig. 7(d). A model resolution analysis determines to what extent the true model maps into an inversion model. The smaller the spread of non-zero entries in the row of the resolution matrix around the diagonal entry is and the higher the diagonal entry is, the better is the model resolution provided by the data (Kalscheuer & Pedersen 2007; Kalscheuer et al.2010, 2013). Due to the limited space, the resolving kernels of only a few cells of the joint inversion model constrained by bathymetry and water resistivities are shown in Fig. 10 and the relevant information of the selected cells as well as the linearized model errors f are listed in Table 2. Underneath the land parts of the profile, the model resistivities are partly resolved to a depth of about 20 m (e.g. cell a in Fig. 10a). However, the main lobe of the resolving kernel is already vertically offset from the investigated cell. The sediments below the lake are also rather well resolved (cell b in Fig. 10b), but below the sediments the corresponding resolving kernels are not focused (cells c, e and d in Figs 10c and d) because the conductive water and sediments together limit the penetration depth considerably. Cells representative of the speculated fracture zones are partly resolved (cells f and g in Figs 10a and b). However, the corresponding resolving kernels show only positive side lobes on the cells under investigation while the main lobes are located at shallow depth (<10 m). Possible reasons for the limited resolution of the speculated fracture zones are noise in the data (the RMT and ERT data are noisy in those parts, see Figs 8 and 9), the general experimental setup (one of our main targets is located close to the NW end of the profile) and fundamental physical equivalences (non-uniqueness of the ill-posed inverse problem). Figure 10. View largeDownload slide Resolving kernels for seven model parameters (a–g) of the joint inversion model constrained by bathymetry and water resistivity measurements in Fig. 7(d). The considered cells are marked by white diamonds in the figure panels. The red crosses depict the centres of resolution and the horizontal and vertical resolution lengths (Kalscheuer & Pedersen 2007). The isolines are the logarithmic resistivity values of the model in Fig. 7(d). Two cells are shown in one subfigure to take less space. Model parts with low resolution show a shift between the investigated cells (diamonds) and the main lobes of the resolving kernels (cell d as an example). Figure 10. View largeDownload slide Resolving kernels for seven model parameters (a–g) of the joint inversion model constrained by bathymetry and water resistivity measurements in Fig. 7(d). The considered cells are marked by white diamonds in the figure panels. The red crosses depict the centres of resolution and the horizontal and vertical resolution lengths (Kalscheuer & Pedersen 2007). The isolines are the logarithmic resistivity values of the model in Fig. 7(d). Two cells are shown in one subfigure to take less space. Model parts with low resolution show a shift between the investigated cells (diamonds) and the main lobes of the resolving kernels (cell d as an example). Table 2. Positions (yc, zc) and extents (Δy, Δz) of the cells chosen for model error and resolution analyses as well as linearized error factors f of the resistivities of these cells. Resolving kernels are depicted in Fig. 10. Cell lable  yc  zc  Δy  Δz  f  a  100.2  13.2  4.9  2.1  1.09  b  199.6  18.5  5.0  3.4  1.10  c  299.4  22.3  5.0  4.3  1.10  d  449.1  27.2  5.0  5.4  1.10  e  499.1  22.3  5.0  4.3  1.10  f  548.7  18.5  5.0  3.4  1.13  g  597.7  13.2  5.0  2.1  1.08  Cell lable  yc  zc  Δy  Δz  f  a  100.2  13.2  4.9  2.1  1.09  b  199.6  18.5  5.0  3.4  1.10  c  299.4  22.3  5.0  4.3  1.10  d  449.1  27.2  5.0  5.4  1.10  e  499.1  22.3  5.0  4.3  1.10  f  548.7  18.5  5.0  3.4  1.13  g  597.7  13.2  5.0  2.1  1.08  View Large 4.8 Geological interpretation The geological interpretation is summarized in Fig. 11(a). Three distinct fracture zones exist at distances of 200–300 m along the profile based on previous studies (Wikberg et al.1991; Stanfors et al.1999). These fracture zones are not well resolved in our models (EW-7 is the only one detected, Fig. 11a) because the overlying conductive sea water and sediments limit the resolution in this part of the profile. Following Wikberg et al. (1991) and Stanfors et al. (1999), the fracture zones NE-3, NE-4 and EW-7 should exist in this part of the profile. In the resistivity models shown, the sediments above the fracture zones are more conductive than the sea water (see also Ronczka et al.2017). At 500–700 m distance along the profile, the model suggests two distinct conductive zones (possible fracture zones). The one from 500–580 m distance is likely related to EW-5 (Wikberg et al.1991), which was previously a poorly constrained fracture zone. However, more information is needed to conclusively confirm the existence of EW-5. The other conductive zone at 600–700 m distance along the profile is most probably related to the NE-1 fracture system. The 65° dip of NE-1 suggested by Berglund et al. (2003) matches closely the models obtained by the joint inversion (Fig. 11a). The position of NE-1 in the inversion models matches the surface projection of NE-1 observed in the tunnel (Figs 3b and 11a). However, the projection could be slightly shifted compared with the real position, because it is projected from the tunnel and limited boreholes to the ground surface and no relevant observation is available on the ground surface. Figure 11. View largeDownload slide (a) Interpretation of joint inversion model in Figure 7(d) with possible positions of fracture zones EW-7, NE-4, NE-3, EW-5 and NE-1. Black thick arrows on top of the model section mark the positions of the fracture zones along the profile based on the previous studies shown in Fig. 3. Black lines at depth mark interpreted boundaries of fracture zones as inferred from the model in panel (a). (b) Test model generated with main structures from joint inversion model constrained by bathymetry and water resistivity measurements at 500–700 m distance and from previous studies at 200–350 m distance (EW-7 could be observed in our model). (c) Single inversion model of synthetic RMT TE-mode data. (d) Single inversion result of synthetic ERT data. (e) Joint inversion model of synthetic RMT TE-mode and ERT data. (f) Table showing RMS. The joint inversion model fits both data sets reasonably well and shows the highest level of detail among all the inversions. Figure 11. View largeDownload slide (a) Interpretation of joint inversion model in Figure 7(d) with possible positions of fracture zones EW-7, NE-4, NE-3, EW-5 and NE-1. Black thick arrows on top of the model section mark the positions of the fracture zones along the profile based on the previous studies shown in Fig. 3. Black lines at depth mark interpreted boundaries of fracture zones as inferred from the model in panel (a). (b) Test model generated with main structures from joint inversion model constrained by bathymetry and water resistivity measurements at 500–700 m distance and from previous studies at 200–350 m distance (EW-7 could be observed in our model). (c) Single inversion model of synthetic RMT TE-mode data. (d) Single inversion result of synthetic ERT data. (e) Joint inversion model of synthetic RMT TE-mode and ERT data. (f) Table showing RMS. The joint inversion model fits both data sets reasonably well and shows the highest level of detail among all the inversions. A synthetic model (Fig. 11b) was designed based on these interpretations to evaluate the suggested resolution properties of the joint inversion model. Synthetic lake-floor ERT and RMT data sets were generated at the same observation positions as the field data and contaminated with noise corresponding to the data error floors used in the inversions of the field data. Furthermore, the same data editing was used in subsequent inversions to better understand the inversions of field data. The results of single inversions (Figs 11c and d) and joint inversion (Fig. 11e) of the synthetic data reflect well how the field data were modeled in the single and joint inversions. The RMT single inversion resolved the facture zone at 600–700 m distance. Because of the low data coverage and structural complexity of the fracture zones at 500–700 m distance, the ERT single inversion has only incorrectly predicted one fracture zone (Fig. 11d). The joint inversion, however, resolves both fracture zones (Fig. 11e). The approximate maximal penetration depth of the combined data sets is also presented in Fig. 11(e) to evaluate the joint inversion results. 5 DISCUSSION AND CONCLUSIONS 5.1 Lake-floor ERT and boat-towed RMT joint inversion and useful techniques For the first time, a joint inversion of boat-towed RMT and lake-floor ERT data has been implemented based on a few modifications of an existing joint inversion code and a well-designed field case study. The implementation was tested with both synthetic and field data sets. All the results demonstrated that the joint inversion of boat-towed RMT TE-mode and lake-floor ERT data performs better than individual inversions. This is mainly due to: (1) lake-floor ERT and boat-towed RMT data together lead to improved data coverage of geological targets. (2) The sensitivities of the two methods complement each other (Candansayar & Tezkan 2008; Kalscheuer et al.2010). (3) The two methods are affected differently by noise. This also explains why a joint inversion model can fit both data sets, but the model inverted from one data set cannot explain the other data set. The inversion code EMILIA offers two different weighting schemes. One allows modifying the weights of different data sets and the other one allows modifying the weights on smoothness constraints in horizontal and vertical directions. Stronger weight on the data set, which has the lower number of data points, combines information from both data sets more evenly and results in a model that can explain both data sets (Fig. 7b). Weighting of smoothness constraints in different model directions is often required to explain data properly and, thus, it may eliminate artefacts produced by inversion with equal weights. In this study, this technique was used to remove possible artefacts in the vertical direction of the model (Fig. 7c). However, this did not influence the conductive zones at 550–700 m distance along the profile, i.e. the fracture system NE-1 and the possible fracture system EW-5 (Fig. 7c), which were already observed in the model in Fig. 7(b). This supports the hypothesis that the fracture zones NE-1 and EW-5 are resolved. A priori information provides known model structure to the inversion, such that ambiguity can be reduced in those parts of the model that are related to the a priori information (Siripunvaraporn et al.2005). This may lead to more accurate estimates of the model parameters that were determined with higher uncertainty without using a priori information. In the joint inversion using a priori information, the water resistivities, obtained from direct measurements, were fixed. Additionally, the smoothness constraints were decoupled along the lower boundary of the water body in order to correctly estimate the resistivity and thickness of the conductive lake sediments. All those setups together improved the inversion results for the fracture systems (Fig. 7d). 5.2 Evaluation of joint inversion and interpretation Exploration depth investigations, data misfit analyses and model resolution and error analyses were our main methods to evaluate the inversion models which had a reasonable RMS. Using only RMS as a standard evaluation tool is far from being sufficient. The well-constrained part of an inversion model is distinguished by being located within the depth of penetration range and represented by data with low misfit, focused model resolution and reasonably low model errors (Figs 7d, 8f and 9f and Tables 1 and 2). In this way, the well-constrained part of the model is deemed to be representative of the subsurface structures. One fracture zone (NE-1) is well delineated in this study (fracture zones EW-5 and EW-7 are also inferred although with a lower level of confidence) using the joint inversion approach applied to the boat-towed RMT and lake-floor ERT data sets. A synthetic test based on the interpretation results and another test just incorporating land RMT data (Fig. 12) proved that the joint inversion of the boat-towed RMT and lake-floor ERT data was indeed capable of resolving fractures in the bedrock. Detailed information about the geometry of three fracture zones (NE-1, EW-5 and EW-7) together with thick sediments at 200–300 m along the profile are new useful findings for SKB to consider in their ongoing research and can be drilled to be confirmed. However, the fracture zones NE-3 and NE-4 below the deep water and thick sediments cannot be resolved, because resolution underneath the conductive water and more conductive sediments is dramatically reduced. The promising results from the case study indicate that using boat-towed RMT and lake-floor ERT together would have higher resolution, if water and sediments were less conductive. Figure 12. View largeDownload slide Comparison between the joint inversion models inverted with boat-towed RMT data and without boat-towed RMT data. (a) Joint inversion with boat-towed RMT data, the figure is also shown in Fig. 6(c). (b) Joint inversion without boat-towed RMT data. (c) Table showing RMS. Both inversions have the same setup with regard to mesh, weighting and smoothness constraints. Apparently, boat-towed RMT data contribute information to the joint inversion model underneath the shallow water parts. Figure 12. View largeDownload slide Comparison between the joint inversion models inverted with boat-towed RMT data and without boat-towed RMT data. (a) Joint inversion with boat-towed RMT data, the figure is also shown in Fig. 6(c). (b) Joint inversion without boat-towed RMT data. (c) Table showing RMS. Both inversions have the same setup with regard to mesh, weighting and smoothness constraints. Apparently, boat-towed RMT data contribute information to the joint inversion model underneath the shallow water parts. Figure B1. View largeDownload slide (a) Inversion model computed from synthetic RMT determinant impedance data. (b) Inversion model computed from joint inversion of synthetic RMT determinant impedance and ERT data. (c) Inversion model computed from synthetic RMT TM-mode data. (d) Joint inversion model for synthetic RMT TM-mode and ERT data. (e) Inversion model for RMT TE- and TM-mode data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. The true model for which the synthetic data were computed is shown in Fig. 2(a). Symbols are the same as in Fig, 2. Figure B1. View largeDownload slide (a) Inversion model computed from synthetic RMT determinant impedance data. (b) Inversion model computed from joint inversion of synthetic RMT determinant impedance and ERT data. (c) Inversion model computed from synthetic RMT TM-mode data. (d) Joint inversion model for synthetic RMT TM-mode and ERT data. (e) Inversion model for RMT TE- and TM-mode data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. The true model for which the synthetic data were computed is shown in Fig. 2(a). Symbols are the same as in Fig, 2. 5.3 Prospect for lake-floor ERT and boat-towed RMT joint inversion With increasing population, higher demands will be imposed on future transportation systems not just on land but also under lakes and rivers. In Sweden, an area of about 40 000 km2 is covered by inland water, constituting around 95 700 lakes. Inevitably, parts of the transportation system have to pass either over or under the water. In the planning phase, site investigations of these areas will benefit from boat-towed RMT to design the optimal profile for lake-floor ERT, because boat-towed RMT is much more efficient than ERT. Then, the boat-towed RMT and lake-floor ERT joint inversion can further be used to map the most interesting area. This joint approach is more capable of resolving weakness zones effectively, for examples, unconsolidated sediments, faults and fracture zones, than single-method investigations based on either ERT or RMT. Thus, joint inversion of lake-floor ERT and boat-towed RMT data is a better way to improve model resolution in investigations of geological structure below shallow water bodies with those two geophysical methods at the present stage. This method can effectively resolve underwater structures and reduce the cost of pre-investigations in underwater infrastructure planning. For deeper (>10 m depending on the case) and conductive (<5 Ωm) water bodies, the RMT method will not provide sufficient depth penetration to explore structures below the water bodies. In such cases, a combination of boat-towed controlled-source RMT (in the frequency range of 1–10 kHz, Wang et al.2017), or boat-towed transient EM with a transmitter loop of adequate size (e.g. Barrett et al.2005; Hatch et al.2010; Mollidor et al.2013; Bekesi et al.2014) and lake-floor ERT methods should be employed. Acknowledgements This work was carried out within the frame of TRUST2.2 (http://trust-geoinfra.se) project sponsored by Formas (project number: 25220121907), SGU, BeFo, SBUF/Skanska, FQM and NGI. The ERT data were acquired within the TRUST 2.1 and 4.2 (http://trust-geoinfra.se) projects funded by Formas, BeFo and SBUF/Skanska. SKB provided access to the site and provided support for the measurements. Nova FoU is thanked for partly funding this work. UPPMAX (project numbers: SNIC 2016-1-38 and SNIC 2016-1-439) is acknowledged for providing computation facilities. M. Ronczka and P. Olsson are thanked for providing the ERT data. A.A. Pfaffhuber, M. Everett and one anonymous reviewer are thanked for their constructive and useful comments. S. Wang thanks the China Scholarship Council (201306370039), Formas (25220121907), National Natural Science Foundation of China (41227803) and National High Technology Research and Development Program of China (2012AA09A20105) for supporting his PhD study at Uppsala University. REFERENCES Abubakar A., Gao G., Habashy T.M., Liu J., 2012. Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data, Inverse Prob. , 28( 5), 055016, doi:10.1088/0266-5611/28/5/055016. Google Scholar CrossRef Search ADS   Almén K.E., Stenberg L., 2005. Characterisation Methods and Instruments . Experiences from the construction phase (no. SKB-TR-05-11). Äspö Hard Rock Laboratory. Swedish Nuclear Fuel and Waste Management Company. Barrett B., Heinson G., Hatch M., Telfer A., 2005. River sediment salt-load detection using a water-borne transient electromagnetic system, J. Appl. Geophys. , 58( 1), 29– 44. Google Scholar CrossRef Search ADS   Bastani M., 2001. EnviroMT—a new controlled source/radio magnetotelluric system, PhD thesis , Uppsala University. Bastani M., Hübert J., Kalscheuer T., Pedersen L.B., Godio A., Bernard J., 2012. 2D joint inversion of RMT and ERT data versus individual 3D inversion of full tensor RMT data: an example form Trecate site in Italy, Geophysics , 77( 4), WB233– WB243. Google Scholar CrossRef Search ADS   Bastani M., Pedersen L.B., 2001. Estimation of magnetotelluric transfer functions from radio transmitters, Geophysics , 66( 4), 1038– 1051. Google Scholar CrossRef Search ADS   Bastani M., Persson L., Mehta S., Malehmir A., 2015. Boat-towed radio-magnetotellurics—a new technique and case study from the city of Stockholm, Geophysics , 80( 6), B193– B202. Google Scholar CrossRef Search ADS   Bekesi G., Telfer A., Woods J., Forward P., Burnell R., Hatch M., 2014. Quantitative measure of salt interception using in-river transient electromagnetic geophysics, Aust. J. Water Res. , 18( 1), 55– 69. Berglund J., Curtis P., Eliasson T., Olsson T., Starzec P., Tullborg E.L., 2003 Update of the Geological Model 2002 . Äspö Hard Rock Laboratory. Swedish Nuclear Fuel and Waste Management Company. Brodic B., Malehmir A., Juhlin C., Dynesius L., Bastani M., Palm H., 2015. Multicomponent broadband digital-based seismic landstreamer for near-surface applications, J. Appl. Geophys. , 123, 227– 241. Google Scholar CrossRef Search ADS   Brodic B., Malehmir A., Juhlin C., 2017. Delineating fracture zones using surface-tunnel-surface seismic data, P-S and S-P mode conversions, J. Geophys. Res. , doi:10.1002/2017JB014304. Candansayar M.E., Tezkan B., 2008. Two-dimensional joint inversion of radiomagnetotelluric and direct current resistivity data, Geophys. Prospect. , 56( 5), 737– 749. Google Scholar CrossRef Search ADS   Christensen N.B., 1998. The determination of electrical anisotropy using surface electric and electromagnetic methods, in 11th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems , doi:org/10.4133/1.2922495. Constable S.C., Parker R.L., Constable C.G., 1987. Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics , 52( 3), 289– 300. Google Scholar CrossRef Search ADS   Coutant O., Bernard M.L., Beauducel F., Nicollin F., Bouin M.P., Roussel S., 2012. Joint inversion of P-wave velocity and density, application to La Soufrière of Guadeloupe hydrothermal system, Geophys. J. Int. , 191( 2), 723– 742. Google Scholar CrossRef Search ADS   Cosma C., Olsson O., Keskinen J., Heikkinen P., 2001. Seismic characterization of fracturing at the Äspö Hard Rock Laboratory, Sweden, from the kilometer scale to the meter scale, Int. J. Rock Mech. Min. Sci. , 38( 6), 859– 865. Google Scholar CrossRef Search ADS   Dahlin T., Zhou B., 2006. Multiple-gradient array measurements for multichannel 2D resistivity imaging, Near Surf. Geophys ., 4( 2), 113– 123. Dahlin T., Loke M.H., Siikanen J., Höök M., 2014. Underwater ERT survey for site investigation for a new line for Stockholm Metro, in Near Surface Geoscience 2014–20th European Meeting of Environmental and Engineering Geophysics , doi:10.3997/2214-4609.20142060. Dey A., Morrison H.F., 1979. Resistivity modelling for arbitrarily shaped two-dimensional structures, Geophys. Prospect. , 27( 1), 106– 136. Google Scholar CrossRef Search ADS   Daily W., Ramirez A., Binley A., LaBrecque D., 2005. Electrical resistance tomography—theory and practice, in Near-Surface Geophysics , pp. 525– 550, ed. Butler D.K., Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Edwards L.S., 1977. A modified pseudosection for resistivity and IP, Geophysics , 42( 5), 1020– 1036. Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near‐surface materials by joint 2D inversion of dc resistivity and seismic data, Geophys. Res. Lett. , 30( 13), doi:10.1029/2003GL017370. Gallardo L.A., Meju M.A., 2004. Joint two‐dimensional DC resistivity and seismic travel time inversion with cross‐gradients constraints, J. Geophys. Res. , 109(B3), doi:10.1029/2003JB002716. Gallardo L.A., Meju M.A., 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification, Geophys. J. Int. , 169( 3), 1261– 1272. Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2011. Structure‐coupled multiphysics imaging in geophysical sciences, Rev. Geophys. , 49( 1), doi:10.1029/2010RG000330. Haber E., Gazit M.H., 2013. Model fusion and joint inversion, Surv. Geophys. , 34( 5), 675– 695. Google Scholar CrossRef Search ADS   Hatch M., Munday T., Heinson G., 2010. A comparative study of in-river geophysical techniques to define variations in riverbed salt load and aid managing river salinization, Geophysics , 75( 4), WA135– WA147. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints, Geophysics , 74( 6), R99– R109. Google Scholar CrossRef Search ADS   Huang H., 2005. Depth of investigation for small broadband electromagnetic sensors, Geophysics , 70( 6), G135– G142. Google Scholar CrossRef Search ADS   Jones A.G., Groom R.W., 1993. Strike-angle determination from the magnetotelluric impedance tensor in the presence of noise and local distortion: rotate at your peril!, Geophys. J. Int. , 113( 2), 524– 534. Google Scholar CrossRef Search ADS   Kalscheuer T., Pedersen L.B., 2007. A non-linear truncated SVD variance and resolution analysis of two-dimensional magnetotelluric models, Geophys. J. Int. , 169( 2), 435– 447. Google Scholar CrossRef Search ADS   Kalscheuer T., Pedersen L.B., Siripunvaraporn W., 2008. Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents, Geophys. J. Int. , 175( 2), 486– 514. Google Scholar CrossRef Search ADS   Kalscheuer T., García Juanatey M.A., Meqbel N., Pedersen L.B., 2010. Non-linear model error and resolution properties from two-dimensional single and joint inversions of direct current resistivity and radiomagnetotelluric data, Geophys. J. Int. , 182( 3), 1174– 1188. Google Scholar CrossRef Search ADS   Kalscheuer T., Bastani M., Donohue S., Persson L., Pfaffhuber A.A., Reiser F., Ren Z.Y., 2013. Delineation of a quick clay zone at Smørgrav, Norway, with electromagnetic methods under geotechnical constraints, J. Appl. Geophys. , 92, 121– 136. Google Scholar CrossRef Search ADS   Kalscheuer T. et al., 2015. Joint inversions of three types of electromagnetic data explicitly constrained by seismic observations: results from the central Okavango Delta, Botswana. Geophys. J. Int. , 202( 3), 1429– 1452. Google Scholar CrossRef Search ADS   Lines L.R., Schultz A.K., Treitel S., 1988. Cooperative inversion of geophysical data, Geophysics , 53( 1), 8– 20. Google Scholar CrossRef Search ADS   Loke M.H., 1999. Electrical imaging surveys for environmental and engineering studies, in A Practical Guide to 2-D and 3-D Surveys . Loke M.H., 2002. Rapid 2D resistivity forward modelling using the finite-difference and finite-element methods. Loke M.H., Lane J.W., 2004. Inversion of data from electrical resistivity imaging surveys in water-covered areas, Explor. Geophys. , 35, 266– 271. Google Scholar CrossRef Search ADS   Makurat A., Løset F., Hagen A.W., Tunbridge L., Kveldsvik V., Grimstad E., 2006. A Descriptive Rock Mechanics Model for the 380–500 m Level . Report R-02-11, SKB. Malehmir A. et al., 2015. Planning of urban underground infrastructure using a broadband seismic landstreamer—tomography results and uncertainty quantifications from a case study in southwestern Sweden, GeoPhysics , 80( 6), B177– B192. Google Scholar CrossRef Search ADS   Mehta S., Bastani M., Malehmir A., Pedersen L.B., 2017. Resolution and sensitivity of boat-towed RMT data to delineate fracture zones–example of the Stockholm bypass multi-lane tunnel, J. Appl. Geophys. , 139, 131– 143. Google Scholar CrossRef Search ADS   Meju M.A., 1996. Joint inversion of TEM and distorted MT soundings: some effective practical considerations, Geophysics , 61( 1), 56– 65. Google Scholar CrossRef Search ADS   Menke W., 1989. Geophysical Data Analysis: Discrete Inverse Theory , International Geophysics Series , Vol. 45. Academic Press, London. Mollidor L., Tezkan B., Bergers R., Löhken J., 2013. Float‐transient electromagnetic method: in‐loop transient electromagnetic measurements on Lake Holzmaar, Germany, Geophys. Prospect. , 61( 5), 1056– 1064. Google Scholar CrossRef Search ADS   Monteiro Santos F.A., Andrade Afonso A.R., Dupis A., 2006. 2D joint inversion of dc and scalar audio-magnetotelluric data in the evaluation of low enthalpy geothermal fields, J. Geophys. Eng. , 4( 1), 53, doi:10.1088/1742-2132/4/1/007. Google Scholar CrossRef Search ADS   Moorkamp M., Heincke B., Jegen M., Roberts A.W., Hobbs R.W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data, Geophys. J. Int. , 184( 1), 477– 493. Google Scholar CrossRef Search ADS   Moorkamp M., Jones A.G., Fishwick S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data, J. Geophys. Res. , 115( B4), B04318, doi:10.1029/2009JB006369. Google Scholar CrossRef Search ADS   Moorkamp M., Lelievre P.G., Linde N., Khan A., 2016. Integrated Imaging of the Earth . Wiley Press. Google Scholar CrossRef Search ADS   Pedersen L.B., Engels M., 2005. Routine 2D inversion of magnetotelluric data using the determinant of the impedance tensor, Geophysics , 70( 2), G33– G41. Google Scholar CrossRef Search ADS   Ronczka M., Hellman K., Günther T., Wisén R., Dahlin T., 2017. Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory, Solid Earth , 8, 671– 682. Google Scholar CrossRef Search ADS   Rhén I., Gustafsson G., Stanfors R., Wikberg P., 1997. Äspö HRL-Geoscientific Evaluation 1997/5. Models Based on Site Characterization 1986–1995  (no. SKB-TR–97-06). Swedish Nuclear Fuel and Waste Management Company. Sasaki Y., 1989. Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data, Geophysics , 54( 2), 254– 262. Google Scholar CrossRef Search ADS   Schlumberger M., 1939. The application of telluric currents to surface prospecting, EOS, Trans. Am. Geophys. Un. , 20( 3), 271– 277. Google Scholar CrossRef Search ADS   Silvester P.P., Ferrari R.L., 1996. Finite Elements for Electrical Engineers . Cambridge University Press. Google Scholar CrossRef Search ADS   Stanfors R., Rhén I., Tullborg E.L., Wikberg P., 1999. Overview of geological and hydrogeological conditions of the Äspö hard rock laboratory site, Appl. Geochem. , 14( 7), 819– 834. Google Scholar CrossRef Search ADS   Sudha A., Tezkan B., Siemon B., 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter‐borne electromagnetic data, Geophys. Prospect. , 62( 3), 597– 614. Google Scholar CrossRef Search ADS   Siripunvaraporn W., Egbert G., 2000. An efficient data-subspace inversion method for 2D magnetotelluric data, Geophysics , 65( 3), 791– 803. Google Scholar CrossRef Search ADS   Siripunvaraporn W., Egbert G., Lenbury Y., Uyeshima M., 2005. Three-dimensional magnetotelluric inversion: data-space method, Phys. Earth Planet. Inter. , 150( 1), 3– 14. Google Scholar CrossRef Search ADS   Spies B.R., 1989. Depth of investigation in electromagnetic sounding methods, Geophysics , 54( 7), 872– 888. Google Scholar CrossRef Search ADS   Swift C.M. Jr, 1967. A magnetotelluric investigation of an electrical conductivity anomaly in the South Western United States, PhD thesis , MIT, Cambridge press. Tezkan B., Goldman M., Greinwald S., Hördt A., Müller I., Neubauer F.M., Zacher G., 1996. A joint application of radiomagnetotellurics and transient electromagnetics to the investigation of a waste deposit in Cologne (Germany), J. Appl. Geophys. , 34( 3), 199– 212. Google Scholar CrossRef Search ADS   Tezkan B., Hördt A., Gobashy M., 2000. Two-dimensional radiomagnetotelluric investigation of industrial and domestic waste sites in Germany, J. Appl. Geophys. , 44( 2), 237– 256. Google Scholar CrossRef Search ADS   TRUST, 2013. Available at: http://trust-geoinfra.se/, (accessed 2016 July 20). Tryggvason A., Linde N., 2006. Local earthquake (LE) tomography with joint inversion for P‐and S‐wave velocities using structural constraints, Geophys. Res. Lett. , 33( 7), L07303, doi:10.1029/2005GL025485. Google Scholar CrossRef Search ADS   Turberg P., Müller I., Flury F., 1994. Hydrogeological investigation of porous environments by radio magnetotelluric-resistivity (RMT-R 12–240 kHz), J. Appl. Geophys. , 31( 1-4), 133– 143. Google Scholar CrossRef Search ADS   Vozoff K., Jupp D.L.B., 1975. Joint inversion of geophysical data, Geophys. J. R. Astr. Soc. , 42, 977– 991. Google Scholar CrossRef Search ADS   Wang S., Bastani M., Kalscheuer T., Malehmir A., Dynesius L., 2017. Controlled source boat-towed radio-magnetotellurics for site investigation at Äspö Hard Rock Laboratory, Southeastern Sweden, in 79th EAGE Conference and Exhibition 2017 , doi:10.3997/2214-4609.201700565. Wikberg P., Gustafson G., Rhén I., Stanfors R., 1991. Evaluation and Conceptual Modelling Based on the Pre-investigations 1986–1990  (no. SKB-TR–91-22). Aespoe Hard Rock Laboratory.  Swedish Nuclear Fuel and Waste Management Company. Xiong B., Wang S.G., 2011. Wave-number domain features of primary field of H− Hz arrangement wide field electromagnetic method, in International Conference on Instrumentation, Measurement, Circuits and Systems (ICIMCS 2011) . ASME Press, doi:10.1115/1.859902.paper196. Xu S.Z., Duan B.C., Zhang D.H., 2000. Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2.5 D electrical modelling, Geophys. Prospect. , 48( 5), 789– 796. Google Scholar CrossRef Search ADS   Yogeshwar P., Tezkan B., Israil M., Candansayar M.E., 2012. Groundwater contamination in the Roorkee area, India: 2D joint inversion of radiomagnetotelluric and direct current resistivity data, J. Appl. Geophys. , 76, 127– 135. Google Scholar CrossRef Search ADS   Zhang P., Roberts R.G., Pedersen L.B., 1987. Magnetotelluric strike rules, Geophysics , 52( 3), 267– 278. Google Scholar CrossRef Search ADS   Zonge K., Wynn J., Urquhart S., 2005. Resistivity, induced polarization, and complex resistivity, in Near-Surface Geophysics , pp. 265– 300, ed. Butler, D.K., Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   APPENDIX A: SENSITIVITY CALCULATION Two types of ERT field data are usually collected in field applications. One is apparent resistivity (eq. 7), which is an average resistivity of the part of the subsurface penetrated by the ERT current system and equal to the real material resistivity when the subsurface is a homogeneous half-space. Another one is resistance (eq. 8). Both of them can be used to invert for a resistivity model. In EMILIA, apparent resistivities are transformed to logarithmic values to stabilize and speed up the inversion. Similarly, model resistivities are transformed to logarithmic values to speed up the inversion and avoid negative resistivity values. Thus, sensitivities are computed for logarithmic apparent resistivities with regard to logarithmic model resistivities yielding the following relationship to the sensitivities of the original quantities:   $$\frac{\partial \ \log _{10} (\rho _{a})}{\partial \ \log _{10} (\rho)} = - \frac{\partial \ \log _{10} (\rho _{a})}{\partial \ \log _{10} (\sigma )} = - \frac{\sigma }{\rho _{a}} \frac{\partial \rho _{a}}{\partial \sigma },$$ (A1)where σ = 1/ρ is conductivity. However, when resistances are used in inversion, it is less meaningful than for apparent resistivities to use logarithmic values, because resistances could be negative. For resistance data, sensitivities are calculated with regard to logarithmic model resistivities establishing the following relationship to the sensitivities of the original quantities   \begin{eqnarray} \frac{{\partial R}}{{\partial \log_{10} (\rho)}} &=&- \frac{{\partial R}}{{\partial \log_{10} (\sigma)}} = - \frac{{\partial R}}{{\partial \left( {\frac{{\ln \sigma }}{{\ln 10}}} \right)}} = - \ln 10\frac{{\partial R}}{{\partial \ln \sigma }}\nonumber\\ &=& - \ln 10 \cdot \sigma \frac{{\partial R}}{{\partial \sigma }}. \end{eqnarray} (A2) EMILIA was modified to use both types of data in inversion. This function makes EMILIA more flexible to utilize field data. APPENDIX B: INVERSIONS OF TM-MODE DATA We have also investigated inversions of synthetic TM mode, TE and TM modes and determinant impedance data computed for the true model in Fig. 2(a). However, the inversion results are only shown in this appendix, because we focus on the results of the TE-mode inversions. In general, the TM mode and determinant impedance single inversions are not in as good structural agreement with the true model as the TE-mode single inversion model, especially the inversion result of the determinant impedances is inferior to that of the TE-mode impedances (Fig. B1a and c). For the joint inversions with the ERT data, good resolution can still be observed in Fig. B1b (joint inversion with RMT determinant impedance data), B1d (joint inversion with RMT TM-mode impedance data) and B1f (joint inversion with RMT TE- and TM-mode impedance data). However, these joint inversion results are not as good as the model from joint inversion of TE-mode and ERT data (Fig. 2d). This is mainly so, because the TM-mode and ERT data have similar directions of electrical current flow. Thus compared to a single inversion of ERT data, the inclusion of TM-mode data may not add significant new information. Since determinant impedance data are geometrical averages of TE- and TM-mode impedances, this effect is also observed in joint inversions of determinant impedance and ERT data to some extent. The joint inversions of the TE- and TM-mode data and of the TE-, TM-mode and ERT data show as good results as the joint inversion of the TE-mode and ERT data. Thus, in our case, the TM-mode data hardly contribute information to our joint inversion of RMT and ERT data (Fig. B1c) that is not already contained in the ERT data. Field data of the TM mode, TE and TM mode, and determinant impedances were also inverted in the case study (TM-mode and tipper data are shown in Fig. B2). In general, the models below the water surface from single inversions of TM mode and determinant impedances (Figs B3a and c) are not as good as that of the TE-mode impedances (Fig. 6a). For the joint inversions, two conductive zones can still be observed in Figs B3(b), (d) and (f). However, they are not as clear as in the model from joint inversion of TE-mode and ERT data (Fig 6c). Figure B2. View largeDownload slide Unedited raw data: (a) TM-mode apparent resistivity data. (b) TM-mode phase data. (c) Real part of vertical magnetic transfer function Ty (Hz/Hy, in our coordinate system). (d) Imaginary part of vertical magnetic transfer function Ty. Figure B2. View largeDownload slide Unedited raw data: (a) TM-mode apparent resistivity data. (b) TM-mode phase data. (c) Real part of vertical magnetic transfer function Ty (Hz/Hy, in our coordinate system). (d) Imaginary part of vertical magnetic transfer function Ty. Figure B3. View largeDownload slide (a) Inversion model computed from RMT determinant impedance field data. (b) Inversion model computed from joint inversion of RMT determinant impedance and ERT field data. (c) Inversion model computed from RMT TM-mode field data. (d) Joint inversion model for RMT TM-mode and ERT field data. (e) Inversion model for RMT TE and TM-mode field data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. For the joint inversion, RMS values of individual data sets are slightly higher than those of the single inversions. Symbols are the same as in Fig. 6. Figure B3. View largeDownload slide (a) Inversion model computed from RMT determinant impedance field data. (b) Inversion model computed from joint inversion of RMT determinant impedance and ERT field data. (c) Inversion model computed from RMT TM-mode field data. (d) Joint inversion model for RMT TM-mode and ERT field data. (e) Inversion model for RMT TE and TM-mode field data. (f) Joint inversion model for RMT TE-mode, RMT TM-mode and ERT field data. (g) Table showing RMS. For the joint inversion, RMS values of individual data sets are slightly higher than those of the single inversions. Symbols are the same as in Fig. 6. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial