Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents

SUMMARY 
 
Electromagnetic surface measurements with the radiomagnetotelluric (RMT) method in the frequency range between 10 and 300 kHz are typically interpreted in the quasi-static approximation, that is, assuming displacement currents are negligible. In this paper, the dielectric effect of displacement currents on RMT responses over resistive subsurface models is studied with a 2-D forward and inverse scheme that can operate both in the quasi-static approximation and including displacement currents. 
 
 
 
Forward computations of simple models exemplify how responses that allow for displacement currents deviate from responses computed in the quasi-static approximation. The differences become most obvious for highly resistive subsurface models of about 3000 Ω m and more and at high frequencies. For such cases, the apparent resistivities and phases of the transverse magnetic (TM) and transverse electric (TE) modes are significantly smaller than in the quasi-static approximation. Along profiles traversing 2-D subsurface models, sign reversals in the real part of the vertical magnetic transfer function (VMT) are often more pronounced than in the quasi-static approximation. On both sides of such sign reversals, the responses computed including displacement currents are larger than typical measurement errors. 
 
 
 
The 2-D inversion of synthetic data computed including displacement currents demonstrates that serious misinterpretations in the form of artefacts in inverse models can be made if displacement currents are neglected during the inversion. Hence, the inclusion of the dielectric effect is a crucial improvement over existing quasi-static 2-D inverse schemes. Synthetic data from a 2-D model with constant dielectric permittivity and a conductive block buried in a highly resistive layer, which in turn is underlain by a conductive layer, are inverted. In the quasi-static inverse model, the depth to the conductive structures is overestimated, artefactual resistors appear on both sides of the conductive block and a spurious conductive layer is imaged at the surface. 
 
 
 
High-frequency RMT field data from Avro, Sweden, are re-interpreted using the newly developed 2-D inversion scheme that includes displacement currents. In contrast to previous quasi-static modelling, the new inverse models have electrical resistivity values comparable to a normal-resistivity borehole log and boundaries between resistive and conductive structures, which correlate with the positions of seismic reflectors.

dielectric effect should be accounted for as soon as the perturbation it causes is roughly equal to the measurement errors of the data. For a typical error level of, say, 2 per cent on the impedance tensor elements, vertically incident plane waves and a relative dielectric permittivity r = 5, it is shown in Section 3.1 that the effect of displacement currents on the impedance phase is above the error level at, for example, f = 15 kHz and ρ = 10 000 m or f = 170 kHz and ρ = 1000 m.
In the following subsections, the existing knowledge of the dielectric effect on plane-wave and controlled-source frequency-domain electromagnetic (FDEM) responses is reviewed. With respect to the RMT method, the plane-wave FDEM responses are of special importance. After that, the resolvability of anomalous dielectric permittivities and previous attempts of quasi-static interpretation of high-frequency RMT data are discussed. In the last part of the introduction, we give an outlook at our 2-D inverse scheme for RMT data that allows for displacement currents and summarize the assumptions we make. Note that for the treatment of the FDEM theory, we choose an exp +iωt time dependence throughout this paper.

Dielectric effect on frequency-domain EM responses
Several publications describe the effect of displacement currents on plane-wave and controlled-source FDEM responses in the VLF and LF frequency ranges, based on analytic solutions by Wait (1953Wait ( , 1970 and Wait & Nabulsi (1996) for a 1-D layered Earth.
In plane-wave FDEM methods like the RMT method, EM fields generated by powerful radio transmitters operating in the VLF and LF frequency ranges are used as primary signals. The aerials employed with the remote radio transmitters are vertical electric dipoles. At distances beyond several free-air wavelengths from the transmitter, that is, in the so-called far-field zone, the EM field essentially resembles that of a plane wave, which is obliquely incident on the Earth's surface (McNeill & Labson 1991). An excellent summary of the theory of plane-wave FDEM impedance, VMT and wave tilt measurements that covers both the quasi-static approximation and the general case with displacement currents, as well as the nature of the radio transmitter source field, is given by Crossley (1981).
For plane-wave EM fields, Sinha (1977) investigates the influence of displacement currents on the wave tilt, that is, the ratio of the horizontal to vertical electric field. On the surface of a homogeneous half-space, both amplitude and phase approach the values of the quasi-static approximation at low frequencies, although they become significantly smaller than the quasi-static responses with increasing frequency.
The dielectric effect on apparent resistivities and phases of radiomagnetotelluric surface impedances is deduced in Crossley (1981), Zacher (1992) and Persson & Pedersen (2002) from 1-D forward computations. On the surface of a homogeneous half-space, both apparent resistivity and phase are smaller than their constant counterparts in the quasi-static approximation. The differences become stronger with increasing frequency. Wait (1953), Sinha (1977), Crossley (1981) and Song et al. (2002) emphasize the importance of the angle of incidence for wave tilt, surface impedance and VMT measurements conducted with plane-wave FDEM methods. The EM field is transmitted vertically into the Earth, independent of the angle of incidence, when the quasi-static approximation is valid. In the general case with displacement currents, however, the angles of incidence and transmission are related through Snell's law. As a consequence, the TM-and TE-mode impedances vary with the angle of incidence and differ at oblique incidence, even if measured on the surface of a layered half-space (Song et al. 2002).
For controlled source air-borne FDEM measurements, Fraser et al. (1990), Huang & Fraser (2002) and Yin & Hodges (2005) simulate responses due to a pair of horizontal coplanar transmitting-receiving coils, operating in the frequency range of 0.4 to 100 kHz. The ratio of secondary magnetic field intensity to primary magnetic field intensity is split into an in-phase component (real part) and a quadrature component (imaginary part). According to Fraser et al. (1990) and Huang & Fraser (2002), displacement currents in the Earth lead to a decrease of the in-phase component and an increase of the quadrature component, compared with the quasi-static case for which both components are positive. The influence of displacement currents in the air (an increase of both components) is rather small compared with that in the Earth (Yin & Hodges 2005).
Using obliquely incident plane waves in the VHF range (30-300 MHz), Nabulsi & Wait (1996) illustrate that a dielectric layer embedded in a highly resistive host is detectable if its thickness and relative permittivity are sufficiently high.
For a controlled source coil-coil FDEM method which operates in the MF (0.3-3 MHz) and HF (3-30 MHz) frequency ranges, Stewart et al. (1994) show that the anomalous response of both a resistive and conductive thin layer is significantly enlarged by the dielectric effect even if there is no contrast of dielectric permittivity between the layers of the model. Stewart et al. (1994) present two field examples, where tilt angle and ellipticity data of the magnetic field polarization ellipse have been successfully inverted for both electric resistivity and dielectric permittivity, with a 1-D inverse scheme.
At frequencies lower than those employed by Nabulsi & Wait (1996) and Stewart et al. (1994), displacement currents become weaker and the resolvability of permittivity anomalies within a limited range of possible relative permittivity values deteriorates. Huang & Fraser (2002) (see Section 1.1) estimate a single value of relative permittivity at their highest frequency of 100 kHz, as it is a badly resolved parameter at lower frequencies. Persson & Pedersen (2002) invert RMT data with frequencies up to 250 kHz for dielectric permittivity, using 1-D models. The differences of inverse models are found to be negligible if the relative dielectric permittivities are limited to the range between 4 and 10, typical of bedrock, and if the resistivities are not larger than 20 000 m (Persson & Pedersen 2002). Relative dielectric permittivities larger than 10 are typical of water bearing sedimentary rocks and soils (Reynolds 1997). Due to the high water content, such formations have relatively low resistivities (typically up to about 500 m), which reduce the importance of displacement currents at VLF and LF frequencies. It is therefore sufficient, in many practical cases of RMT data interpretation, to account for displacement currents by selecting a dielectric permittivity representative of high-resistivity structures in the subsurface.

Quasi-static interpretation of high-frequency RMT data
The difficulties of the interpretation of high-frequency RMT field data in the quasi-static approximation are discussed by Persson & Pedersen (2002) and Linde & Pedersen (2004). For synthetic 1-D RMT impedance responses computed with displacement currents, Persson & Pedersen (2002) compare 1-D inversion results from inverse schemes that utilize both the quasi-static approximation and displacement currents. For a homogeneous half-space model, neglecting displacement currents during the inversion leads to an inverse model with a conductor close to the surface, followed by alternating layers of high and low resistivity at depth (Persson & Pedersen 2002). Similarly, Linde & Pedersen (2004) observe for quasi-static 1-D inversions of RMT field data from the islandÄvrö, Sweden, that a conductive surface layer is modelled more conductive and the underlying unfractured bedrock is modelled more resistive than in the 1-D inversions with displacement currents. The models, due to inversion with displacement currents, are supported by logging data of Gentzschein et al. (1987).
In fact, the work presented by Linde & Pedersen (2004) is a typical example of the interpretation strategies chosen until now, in cases where the dielectric effect in RMT data is to be accounted for. In the absence of a 2-D inversion program that allows for displacement currents, the data interpretation has, so far, been restricted to 1-D inversions with modified analytic forward and Frechet derivative routines, the exclusion of the higher frequency data in 2-D inversions and 3-D forward modelling with the integral equation code X3D by Avdeev et al. (2002).

2-D inversion of RMT data allowing for displacement currents
For the first time, we take displacement currents in a 2-D forward and inverse modelling scheme for RMT data into account by selecting a value of dielectric permittivity that is typical of the subsurface and assuming vertically incident plane waves. As the EM field from remote VLF transmitters can be expected to be incident at an angle closer to 90 • (grazing incidence), it is shown in Section 2.2 that the presence of a moderately resistive surface layer reduces the influence of the angle of incidence considerably. We investigate the effect of displacement currents on 2-D forward responses in the TM-mode, the TE-mode and the VMT and compare our results with the responses computed by the integral equation code X3D by Avdeev et al. (2002), which, at the time of writing, was the only forward code known to us that operates in two or three dimensions and includes displacement currents. Especially, the effect on VMT responses was not considered in the past (cf. Avdeev et al. 2002;Persson & Pedersen 2002). Possible misinterpretations, in the form of artefacts with excessively extreme resistivities in models from quasi-static inverse schemes, are highlighted. The RMT data fromÄvrö (Linde & Pedersen 2004) are re-interpreted with the inverse scheme that allows for displacement currents. The resulting inverse models are compared with the borehole data of Gentzschein et al. (1987) and the seismic reflection model of Juhlin & Palm (1999).
We have added our forward and sensitivity routines, which allow for displacement currents, to the popular 2-D magnetotelluric inverse code REBOCC by Siripunvaraporn & Egbert (2000).

Electromagnetic equations
Assuming a volume of conductivity σ , dielectric permittivity and vacuum permeability μ 0 , Maxwell's equations are written in the frequency domain as displacement current density iω E describes the dielectric effect due to electronic, atomic, molecular and space charge derived polarization of matter with dielectric permittivity in the presence of a time-varying electric field (Keller 1987). In the case that conduction currents dominate over displacement currents (i.e. σ ω ), displacement currents may be neglected in eq. (2). This simplification is known as the quasi-static approximation.
In the following, it is assumed that plane waves are obliquely incident on the Earth's surface in the y-z plane and that the x-direction is the geoelectrical strike direction. Therefore, the admittivityŷ and the EM field components vary only in y and z direction. This choice leads to the definition of the transverse electric (TE) and transverse magnetic (TM) modes for which the vertical electrical and vertical magnetic field components, respectively, vanish. The sets of equations for the TE-and TM-modes are (1) TE-mode: (2) TM-mode: An illustration of the EM field components of the TM-mode and a 2-D subsurface with a cylindrical structure of anomalous electrical properties and infinite extension along the x-axis, that is, the strike direction, is given in Fig. 1. The EM field is obliquely incident at an angle θ 0 , thereby having a wavenumber vector k 0 = (0, k 0y , k 0z ). According to the definition of the TM-mode, the incident, reflected and transmitted magnetic fields H i = (H i x , 0, 0), H r = (H r x , 0, 0), and H t = (H tx , 0, 0), respectively, are all directed along the strike direction whereas the incident, reflected and transmitted electric fields E i = (0, E iy , E iz ), E r = (0, E r y , E rz ), and E t = (0, E ty , E tz ), respectively, are all directed perpendicularly to the strike direction.
In the quasi-static approximation of the TM-mode, j = (0, σ E y , σ E z ) vanishes in the air half-space (Brewitt-Taylor & Weaver 1976) where σ air = 0 is assumed. As a consequence of eqs (9) and (10), H x is then constant in the air half-space, and an inclusion of the air half-space in the modelling domain can be omitted. If displacement currents are accounted for, the magnetic field in the air is no longer independent of the resistivity distribution in the Earth, as the vertical component of the current density is continuous at the air- n z x y Figure 1. EM field components of the TM-mode on a 2-D earth model. The model consists of a structure with anomalous electrical properties that has its strike direction parallel to the x-axis. The EM field is obliquely incident at an angle θ 0 , thereby having a wavenumber vector k 0 = (0, k 0y , k 0z ) = (0, k 0 sin θ 0 , k 0 cos θ 0 ). The incident, reflected and transmitted magnetic fields H i = (H ix , 0, 0), H r = (H rx , 0, 0) and H t = (H tx , 0, 0), respectively, are all directed along the strike direction. The incident, reflected and transmitted electric fields E i = (0, E iy , E iz ), E r = (0, E ry , E rz ), and E t = (0, E ty , E tz ), respectively, are all directed perpendicular to the strike direction. On top of a conductive subsurface, the electromagnetic field is refracted towards the normal, that is, θ t < θ 0 . and in the air, differ from zero (cmp. eqs 9 and 10). Hence, the air half-space must be included in the simulation of the TM-mode. The electric and magnetic field components E x and H x of the TE-and TM-modes, respectively, fulfil the Helmholtz equations (cf. eqs 1 and 2) and where the 2-D assumption ∂/∂x = 0 and eq. (4) were used. In a homogeneous volume, for instance, the general solution of the scalar Helmholtz equations (eqs 13 and 14) is given by Here, k y and k z are the horizontal and vertical components of the wavenumber vector k (see above). The substitution of eq. (15) into eq. (13) yields where k y = k sin θ and k z = k cos θ . The complex wavenumber can be split as k = α − iβ where the real numbers α and β represent propagation and attenuation, respectively, and The inverse of the imaginary part gives the skin depth δ = 1 β over which the amplitude of the EM field is reduced by a factor 1/e. In the quasi-static approximation, the real and imaginary parts are equal, that is, α = β = ωμ 0 σ 2 . The reflection and refraction of plane EM waves at the Earth's surface are governed by Snell's law and the Fresnel equations (Ward & Hohmann 1987). Hence, the EM field measured on the Earth's surface depends on the angle of incidence (see Fig. 1). Three cases of the angle of incidence θ 0 are distinguished. The cases θ 0 = 0 • and θ 0 = ±90 • are known as normal (or vertical) incidence and grazing (or parallel) incidence, respectively. The cases 90 • > θ 0 > 0 • and 0 • > θ 0 > −90 • are called oblique incidence. The refraction of obliquely incident EM waves into the subsurface is conveniently demonstrated for a layered half-space. As a consequence of the boundary conditions for the EM field components at layer interfaces, the horizontal component of the wavenumber vector is constant (Ward & Hohmann 1987), that is, Here, k 0 = ω 2 μ 0 0 is the wavenumber of the air and θ 0 is the angle of incidence. Similarly, k j and θ j are the wavenumber and angle of transmission of the jth layer, respectively. According to eqs (16) and (19) the vertical wavenumber of the jth layer has the form At sufficiently low frequencies, that is, when the quasi-static approximation is valid, k 2 0 /k 2 j = ω 2 μ 0 0 /(ω 2 μ 0 j − iωμ 0 σ j ) → 0 as ω → 0 and k z, j ∼ = k j . Hence, it is only in the quasi-static approximation or at vertical incidence that the EM field is transmitted vertically into the Earth. At high frequencies and oblique incidence, the angle of transmission generally deviates from 0 • .
After solving the Helmholtz equations (eqs 13 and 14) for E x or H x of a 2-D conductivity distribution, the auxiliary fields H y and H z or E y and E z can be computed with eqs (6) and (7) or eqs (9) and (10), respectively.
The off-diagonal elements of the complex 2-D impedance tensor relate the horizontal magnetic fields to the horizontal electric fields of the TE-and TM-mode as and yield the responses commonly used in radiomagnetotellurics, that is, the apparent resistivities and phases φ xy = arg(Z xy ) and φ yx = arg Z yx . Ward et al. (1968) establish a more direct link to the electrical properties of the subsurface, in the general case with displacement currents, by defining an apparent conductivity and an apparent dielectric permittivity. Due to the dependence of the EM field on the angle of incidence, the amplitude and phase of the impedances of the TM-and TE-mode differ even if measured on the surface of a layered Earth. Only if the quasi-static approximation is valid or if the EM field is vertically incident, the TE-and TM-mode impedances of a layered half-space satisfy the relationship Z xy = −Z yx .
For plane waves vertically incident on the surface of a homogeneous half-space with impedivityẑ and admittivityŷ, the TM-mode impedance has the form Z yx = ẑ/ŷ (Wait 1970;Ward & Hohmann 1987). In the quasi-static approximation, the latter expression simplifies to Z yx = √ iωμ 0 /σ , and only in this case, the apparent resistivities and phases measured on a homogeneous half-space equal the resistivity of the half-space and 45 • , respectively.
In the TE-mode, the vertical magnetic field H z is related to the horizontal magnetic field H y through the complex 2-D VMT B: For plane waves obliquely incident on a layered Earth, the VMT generally differs from zero. However, for vertically incident plane waves or in the quasi-static approximation, a VMT that differs from zero is only observed if the admittivityŷ varies laterally (see eqs 5-7).

Normal and oblique incidence
In the case of grazing or oblique incidence, both the incident electric and the incident magnetic fields can have vertical components (see Fig. 1). Already for a 1-D earth model, the TE-and TM-mode are then defined, by demanding that either the electric or the magnetic field be perpendicular to the plane of incidence (Wait 1970;Ward & Hohmann 1987), and the impedance tensor and VMT measured on the Earth's surface depend on the angle of incidence (see Section 2.1). It is therefore important to appraise the error made by assuming vertical incidence during the modelling process. For a layered earth model, the deviations of the TE-and TM-mode impedance amplitudes and phases at an arbitrary angle of incidence from those at normal incidence can be estimated with well-known recurrence formulae (see e.g. Wait 1953Wait , 1970Crossley 1981;Ward & Hohmann 1987;Song et al. 2002).
For the half-space model shown in Fig. 2(a), consisting of two layers with resistivities of 600 and 30 000 m and layer thicknesses of 25 and 75 m, a confining half-space with a resistivity of 600 m and a constant relative permittivity r = 6, the deviations of the amplitude and phase of the TM-and TE-mode impedances from their respective values at normal incidence are shown in Fig. 2. The maximal deviations of 1.5 per cent and 1 • for the amplitude and phase, respectively, occurring at parallel incidence, are of the order of typically expected error levels. A similar model that consists of the uppermost layer underlain by a confining half-space of 30 000 m shows maximal deviations of 1.0 per cent and 0.25 • , respectively, indicating that a considerable part of the distortion in the first case is due to the reflection of the EM energy on the top of the confining half-space.
The angle of incidence can be estimated with the scheme by Song et al. (2002), which requires that the horizontal EM field components are measured simultaneously at adjacent receiver sites. In a typical RMT field campaign, however, a single receiver is moved along the profile. The interpretation is further complicated, as the EM fields of different transmitters, with frequencies close to a nominal frequency, are used to estimate the TM-and TE-mode impedances (Bastani & Pedersen 2001). Generally, the transmitters are off the profile or strike direction and have different angles of incidence; but the angle of incidence, normally, is close to 90 • (grazing incidence) at the site of investigation (Crossley 1981).
As the aerials employed by the remote radio transmitters, typically, are vertical electric dipoles, the incident EM field is that of a TM-mode. Hence, the definitions of TE-and TM-mode based on the geoelectrical structure of the subsurface and on the nature of the incident field are conciliable only for the TM-mode, given that the direction to the remote radio transmitter coincides with the profile direction (as in Fig. 1 for instance). If the transmitter was located off the profile direction, the wavenumber vector k would have an x-component, which, in the general case, would persist within the Earth and invalidate the 2-D assumption ∂/∂x = 0. However, even this problem is amended if a Relative amplitude and phase difference for the TM-and TE-mode impedances with respect to the case of normal incidence for angles of incidence between 0 • (normal incidence) and 90 • (grazing incidence) and at frequencies between 10 and 300 kHz (panels b-e). The earth is assumed to consist of two layers with resistivities of 600 and 30 000 m and layer thicknesses of 25 and 75 m, respectively, a confining half-space with a resistivity of 600 m and a constant relative permittivity r = 6 (panel a). The deviations from the impedance values at normal incidence are largest at grazing incidence.
moderately resistive or conductive surface layer is present, as the EM field is then transmitted almost vertically into the subsurface, and the definition of different modes can be based on the geoelectrical structure. We consider only vertically incident plane-wave fields. As the above example shows, the presence of a moderately resistive or conductive near-surface layer reduces the importance of the angle of incidence, and deviations of the responses for different angles of incidence are then rather small.

Computation of forward responses and sensitivities
The forward problem, that is, the computation of responses for a given model, is solved by discretizing the modelling domain with the finite-difference approximation (FDA), following Hohmann (1987) and Aprea et al. (1997). The derivations of the FDAs for the TE-and TM-modes can be found in Appendix A. Both direct and iterative solvers for the system of linear equations, arising from the FDA of the forward problem, are discussed in Appendix B. As we have not yet managed to implement an appropriate iterative solver, we rely on the LU-decomposition (also known as Gaussian elimination) by Anderson et al. (1999).
The sensitivity matrix J ∈ R N * M describes the perturbations ensuing for N forward responses F[m] ∈ R N due to perturbations of M model parameters m ∈ R M . The entry of the sensitivity matrix for the kth datum with respect to the lth model parameter is then calculated as a partial derivative: The entries of the sensitivity matrix are typically given for the logarithms of the apparent resistivities and the phases of the impedance tensor elements and the real and imaginary parts of the VMT with respect to the logarithms of the cell resistivities. The logarithms are typically chosen relative to the base 10. The sensitivity matrix is computed with the scheme by Rodi (1976) and depends on the FDA of the forward problem. Further information on this algorithm is given in Appendix C. An example of sensitivity matrix entries is given at the end of Section 3.2.

Mesh design
To obtain accurate modelling results, the total extent of the modelling domain (i.e. the finite-difference mesh) and the sizes of individual cells of the finite-difference mesh need to be well adapted to the settings of the experiment, that is, the length of the profile on which measurements were conducted, the lowest and highest frequencies of the measurements and the distributions of electrical conductivity and dielectric permittivity present in the model.
The horizontal and, below the air-Earth interface, the vertical extents of the finite difference mesh must be larger than those used in the quasi-static approximation, as the skin depth δ = 1 β computed with displacement currents (see eq. 18) is larger than its quasi-static counterpart.
Furthermore, the node spacing must be small compared with the scale lengths across which the EM fields vary, that is, the inverse real and imaginary parts of the complex wavenumber k. In the quasi-static approximation, this leads to the well-known requirement that the node spacing must be small compared with the local skin depth (Aprea et al. 1997). In the general case, 1/α < 1/β and the local node spacing must be considerably smaller than 1/α.
A small vertical node spacing is essentially important for the air half-space since the vertically incident plane wavefield propagates undamped (assuming σ air = 0 S m −1 ). In the air, the largest vertical mesh cell dimension must be smaller than 1/α of the highest frequency. This results in the following comparison. In the REBOCC inverse scheme (Siripunvaraporn & Egbert 2000), the conductivity of the air half-space is assumed to be σ air = 10 −10 S m −1 , and the quasi-static skin depth at a frequency of 300 kHz is 92 km. In the general case with displacement currents, σ air = 0 S m −1 and the inverse real part of the wavenumber is 1/α = 159 m for f = 300 kHz. In the former case, the  (2) and (6) in panel i] to the sides of the conductive block. The corresponding maximum and minimum are marked by labels (1) and (7), respectively, in panel (i). The 2-D FDA and integral equation solutions are in good agreement (right-hand column). skin depth exceeds the size of the modelling domain by far, and it is therefore appropriate to address the primary field of the quasi-static case as a uniform inducing field rather than a plane wave incident on the Earth's surface.

Forward modelling examples
We consider a forward modelling example for a homogeneous half-space, with a resistivity of 10 000 m and a relative permittivity r = 5. Assuming vertically incident plane waves, analytic 1-D solutions with the algorithm by Persson & Pedersen (2002) and 2-D FDA solutions were computed for the apparent resistivities and phases of the TM- (Fig. 3) and TE-mode (not shown) at frequencies between 10 and 250 kHz. The comparison of the analytic 1-D solution (marked by a solid line) and the 2-D FDA solution (marked by crosses) shows excellent agreement. At high frequencies, the effect of displacement currents is to decrease the apparent resistivity and phase below the apparent resistivity of 10 000 m and phase of 45 • , respectively, typical of the quasi-static approximation. For a typical error level of 2 per cent on the impedance, the deviations from the quasi-static values are as large as the given errors at 105 kHz for the apparent resistivity and 15 kHz for the phase.
For the simple 2-D model with a block of ρ = 1000 m in a half-space, with a resistivity of ρ = 10 000 m and r = 5 throughout, shown in Fig. 4, 2-D FDA forward responses with displacement currents are compared with both 2-D FDA forward responses for the quasi-static approximation and the 3-D integral equation solution by Avdeev et al. (2002). Responses were computed for the TM-mode impedance, the TE-mode impedance and the VMT. Fig. 5 shows the 2-D FDA forward responses, computed with and without displacement currents in the left-hand column, and the comparison of 2-D FDA forward responses, computed with displacement currents, and 3-D integral equation solutions, with displacement currents, in the right-hand column. The latter comparison indicates that the finite-difference forward scheme is rather accurate. For the given mesh discretization, the relative deviations between the impedance responses of the FDA and integral equation solutions are below 3.0 per cent. The absolute deviations between the VMT responses of the FDA and integral equation solutions are below 0.003. As errors in the computation of two field components might cancel when taking their ratio, a further comparison was done for the 2-D FDA and 3-D integral equation solutions of individual field components (not shown). After an appropriate normalization, the scaled complex field components of the TE-mode deviate by less than 0.7 per cent, whereas the field components of the TM-mode differ by as much as 3.0 per cent. As we do not have insight into the code by Avdeev et al. (2002), it is difficult to give an explanation for the discordance in the latter case.
For  (2) and (6) mark two such pairs of zero transitions in the real part of the VMT and maxima of the imaginary part of the impedance. and 250 kHz (solid lines and star symbols), the influence of displacement currents is considerable, given the chosen resistivity distribution and relative dielectric permittivity.
For stations located on the sides of the conductive block, the effect of displacement currents on TE-and TM-mode impedances is most obvious. Towards the left-and right-hand edges of the mesh, the apparent resistivities and phases approach those of the corresponding homogeneous half-space (see Fig. 3). Also at sites above the conductive block, apparent resistivity and phase are generally smaller than in the quasi-static approximation.
An important effect of displacement currents on the real and imaginary parts of the VMT at high frequencies is the occurrence of lateral sign reversals, located symmetrically around the conductive block. For f = 250 kHz, lateral sign reversals are shown at 260 and 540 m along the profile in the real part of the VMT [marked by labels (2) and (6), respectively, on the solid line in Fig. 5i] and at 75 and 725 m along . Real (upper row) and imaginary (lower row) parts of the normalized current density j x /E x0 at f = 250 kHz of the TE-mode for the general case with displacement currents (left-hand column) and the quasi-static case (right-hand column) for the model shown in Fig. 4. Normalized current densities, which are close to 0 A V −1 m −1 , are plotted in white. Different colourscales were used for the real and imaginary parts. In the general case with displacement currents, the current system penetrates deeper into the subsurface than in the quasi-static case. the profile in the imaginary part of the VMT (solid line in Fig. 5k). In addition to the lateral sign reversals, the real part of the VMT has a maximum at y = 180 m [marked by label (1) in Fig. 5i] and a minimum at y = 620 m [marked by label (7) in Fig. 5i]. The responses at the maximum and minimum are | e(B)| = 0.05. Sign reversals to the sides of the conductive block can also be observed in the real part of the quasi-static response at y = 150 m and y = 650 m (star symbols in Fig. 5i). However, the quasi-static response is comparatively small at sites further away from the block (no larger than | e(B)| = 0.006) and would most likely be masked by noise effects (a typical absolute error is e.g. e(B) ≈ 0.01) if measured in the field. In the general case with displacement currents, the deduction of the horizontal centre of conductive structures from the positions of zero transitions of the VMT B becomes intricate in more complex geological settings. Artefacts might be introduced to inverse models in a quasi-static interpretation.
It is instructive to relate the lateral sign reversals of the VMT to the gradient of the TE-mode impedance Z xy by considering eqs (7) and (21): which corresponds to the following relationships for the real and imaginary parts of the VMT In the synthetic example for f = 250 kHz, the variation of H y away from its approximate 1-D values at the beginning and end of the profile is less than 22 per cent in the quasi-static case (not shown) and less than 12 per cent in the general case (not shown). The real and imaginary parts of the VMT are in good agreement with the expected variation with the lateral derivative of the TE-mode impedance Z xy for the quasi-static (Figs 6c and d) and general cases (Figs 6a and b). For the general case, the positions of the lateral sign reversals in the real part of the VMT and the corresponding maxima in the imaginary part of the impedance are marked with the labels (2) and (6)   (1) (2) (3) (4) (5) (6) (7)  (1) and (7), the magnetic field H h due to currents in the resistive host is larger than the magnetic field H b due to currents in the conductive block, leading to a maximum and a minimum of the real part of the VMT at (1) and (7), respectively. In the vicinity of the conductive block, the magnetic field is dominated by H b resulting in a minimum and a maximum of the real part of the VMT at (3) and (5), respectively. At positions (2) and (6), the vertical components of H h and H b are equal in magnitude but opposite in direction and, hence, the VMT B is zero. The lateral position of the sign reversal at (4) coincides with the centre of the conductive block.
A more quantitative explanation for the lateral zero transitions can be arrived at by investigating eq. (5). As the curl operator treats the real and imaginary parts of H separately, are directly related to the real and imaginary parts of the current density j x =ŷ E x of the TE-mode. However, the current densities of the quasi-static and general cases in the subsurface are not directly comparable. As the propagation of the electric field in the air is modelled differently (i.e. through conduction currents in the quasi-static approximation with a conductivity σ air = 10 −10 Sm −1 and through displacement currents in the general case), there is a large difference in the scale lengths over which the electric field varies in the air (see Section 2.4). This leads to different phases and amplitudes of the electric fields of the two cases at the air-Earth interface, even if equal amplitudes and phases of the electric field are chosen as boundary conditions on the upper edge of the finite-difference mesh. In addition, different vertical node spacings were chosen in the air half-space for the quasi-static and general cases, according to the considerations in Section 2.4. To circumvent this problem, the electric field is scaled by its surface value at the left-hand edge of the mesh. For the general case, the real and imaginary parts of the normalized current density at f = 250 kHz are shown in Figs 7(a) and (c), respectively. Similarly, for the quasi-static case, the real and imaginary parts of the normalized current density at f = 250 kHz are shown in Figs 7(b) and (d), respectively. The area of the highest normalized current density amplitude (up to 5.7 × 10 −4 A V −1 m −1 ) coincides with the conductive block. To the sides of the block at y < 340 m and y > 460 m, the normalized current density amplitude reaches 1.8 × 10 −4 A V −1 m −1 , with only small lateral changes of the real and imaginary parts at the beginning and end of the profile.
An important simplification ensues for the real part of the VMT of the general case, as e(H y ) exceeds m(H y ) by at least a factor 4.4 at all positions along the profile (not shown). Hence, the real part of the VMT can be approximated as e(B) ≈ e(H z )/ e(H y ) and is mostly determined by e( j x ) in Fig. 7(a). We illustrate the sign reversals in the real part of the VMT for the general case with a sketch (Fig. 8a) that describes the real part of the current system in the subsurface and the emerging magnetic field. The real part of the magnetic field due to currents in the resistive host is designated as H h and the part due to currents in the conductive block is designated as H b (Fig. 8a).
To facilitate a simpler comparison, the real part of the VMT for the general case (solid line in Fig. 5i) is plotted in Fig. 8(b). At the beginning and end of the profile [i.e. to the left-hand side of position (1) and to the right-hand side of position (7) in Fig. 8], the lateral homogeneity of the current system generates a magnetic field with a very small H z -component. At positions (1) and (7), that is, at y = 180 and 620 m, the magnetic field H h due to the resistive host is larger than the magnetic field H b due to the conductive block, leading to a maximum and a minimum of the real part of the VMT at (1) and (7), respectively. At positions (2) and (6) to the sides of the block, that is, at y = 260 and 540 m, the H z -components of H h and H b are equal in amplitude but point in opposite directions, leading to zero-transitions of the real part of the VMT. The minimum, zero transition and maximum of the real VMT response at positions (3), (4) and (5), respectively, are similar in both the quasi-static and general cases (see Fig. 5i). Though in magnitude smaller than the current system in the block, the lateral current system is strong enough to generate a commensurable maximum and minimum of the real part of the VMT at y = 180 and 620 m, respectively (Fig. 8b). Hence, the main effect of displacement currents on the real part of the VMT is to increase the response at the edges of the conductor. As noted before, there are no such distinct maxima or minima associated with the lateral sign reversals in the quasi-static VMT response at f = 250 kHz (star symbols in Fig. 5i). The reason is most likely that the vertical extent of the current systems and the total current strengths to the sides of the block (Figs 7b and d) are smaller than in the general case with displacement currents (Figs 7a and c), whereas the current within the conductive block has a comparable amplitude in both cases. It should also be noted that the imaginary part of the VMT increases quite strongly in amplitude if displacement currents are included (Fig. 5k). An explanation with regard to the imaginary part of the current density (shown in Fig. 7c) does not appear to be possible as the imaginary part of H z and the real part of H y are involved.

Inverse modelling examples
Synthetic responses of a simple 2-D model (Fig. 9), with constant relative dielectric permittivity r = 5 and an elongated block with a resistivity of 1000 m that is buried in a resistive layer of 10 000 m and underlain by a half-space of 500 m, were computed for the TMand the TE-mode. The responses were computed at 20 receiver sites for 15 frequencies, ranging from 10 to 250 kHz giving a total of 600 data points. Gaussian white noise, corresponding to 2.5 per cent of the modulus of the computed impedances, was added to the forward responses of both polarizations. After that, two inversions of the synthetic data set were performed with the REBOCC inverse scheme (Siripunvaraporn & Egbert 2000). During the first inversion, displacement currents were allowed for, whereas they were neglected during the second inversion. In both inversions, the error floor was assumed to correspond to 2.5 per cent of the modulus of the impedances, and the starting model was a homogeneous half-space of 10 000 m.
After six iterations with the inversion that allows for displacement currents, a model was obtained (Fig. 10), which fits the data to a rms misfit of 1.04. Additional iterations with REBOCC did not decrease the rms misfit further. The inverse model reproduces the edges of the block and the resistivities of the block and layered half-space rather accurately. The transition from the lower edge of the conductor into the resistive layer is, however, smeared out. Neglecting displacement currents results in convergence problems and an inverse model with many artefactual structures (Fig. 11). The lowest rms misfit of 1.95 was obtained after four iterations. Clearly, an artefactual thin conductive layer is visible at the surface (a similar conductive layer is also observed by Persson & Pedersen (2002) in 1-D inverse models, computed in the quasi-static approximation, for synthetic data of a homogeneous half-space). The lateral extent of the conductive block and the top of the central parts of the block are grossly in error. Two artefactual resistors with resistivities close to 100 000 m appear to the left-and right-hand side of the block. The depth to the top of the underlying conductive layer is shifted from 105 to 130 m. If the synthetic data were generated from a model without the underlying conductive layer, the artefactual resistors would be observed, both to the sides of and below the conductive block (not shown).
A comparison of the relative errors, that is, the differences between the synthetic data and the forward responses scaled by the data errors, generated by the two inverse schemes, is shown in Fig. 12. The relative errors from the inversion that accounts for displacement currents Figure 12. Relative errors obtained from the inversions of the synthetic data of the model shown in Fig. 9. The relative errors from inversions that account for displacement currents and that neglect displacement currents are shown in the left-and right-hand columns, respectively. The relative errors of log 10 ρ a and φ for the TM-mode are shown in panels (a)-(d), whereas the relative errors of log 10 ρ a and φ for the TE-mode are shown in panels (e)-(h). Systematic deviations from the synthetic data are mostly observed at high frequencies and stations to the sides of the conductive block for the model from the inverse scheme that does not allow for displacement currents (Fig. 11) Figure 13. Sensitivity matrix entries for the block model in Fig. 9 and the TE-mode apparent resistivity at f = 250 kHz and a receiver site at y = 90 m (to the left-hand side of the conductive block). The receiver site is marked by a triangle. The edges of the conductive block and the interface between the upper layer and the confining half-space are depicted as solid black lines. Sensitivity values, which are close to 0, are plotted in white. As expected, sensitivities are largest at the left-hand upper edge of the conductive block and the sensitivity entries of the general case (panel a) encompass a larger volume with non-zero values (observe the negative sensitivity values marked in green) than those calculated in the quasi-static approximation (panel b).
(left-hand column of Fig. 12) show relatively random deviations of the forward data from the synthetic data. In contrast to this, the relative errors from the quasi-static inversion (right-hand column of Fig. 12) exhibit systematic deviations in the form of frequency ranges common to groups of neighbouring stations, with relative errors that have absolute values significantly larger than one and the same sign. The systematic deviations originate from the false assumption that displacement currents can be neglected during the inversion. As expected, the misfit is most severe at high frequencies and receiver sites to the sides of the conductive block.
As an example, the row of the sensitivity matrix for the block model in Fig. 9 and the TE-mode apparent resistivity at f = 250 kHz and a receiver site at y = 90 m (to the left-hand side of the conductive block) is shown in Fig. 13. Model parameters with sensitivities close to zero (shown in white colours in Fig. 13) have little influence on the considered data item. As expected, sensitivities, which were computed for the general case (Fig. 13a), encompass a larger volume with non-zero values than those computed in the quasi-static approximation (Fig. 13b). Especially, the depth extend for the non-zero sensitivity values of the general case is larger. This larger depth range is equivalent to a larger depth of investigation for the general case as already indicated in Section 2.4. Linde & Pedersen (2004) investigate highly resistive granitic bedrock on the small islandÄvrö, Sweden, with tensor RMT, in the frequency range of 14-226 kHz. RMT data were acquired on an east-west profile, with a total length of 960 m and a station spacing of 10 m. OnÄvrö, the typical soil thickness is between 0 and 1 m. The bedrock consists mostly of granite. In some locations, aplitic and pegmatitic dykes are encountered (Gentzschein et al. 1987). Previous geophysical studies include borehole measurements by Gentzschein et al. (1987) and a seismic reflection study on the same profile by Juhlin & Palm (1999). A normal-resistivity log and a fracture frequency log of borehole KAV01, located in the central part of the profile, reveal an upper weathered layer, with a thickness of up to 30 m and a resistivity of about 600 m, followed by almost intact and highly resistive bedrock down to a depth of 200 m and with a resistivity between 32 000 and 40 000 m (Gentzschein et al. 1987). Between 200 and 400 m depth, the resistivity slowly decreases to 10 000 m. At greater depth, the bedrock is more fractured and saline pore fluids decrease the electrical resistivity to a few thousand m. Juhlin & Palm (1999) describe two major seismic reflectors (see Fig. 14d) for the depth range down to 400 m. Reflector C is located beneath the western part of the profile, at a depth between 100 and 320 m and dips approximately 60 • to the east. Reflector D is located beneath the central part of the profile at a depth between 150 and 200 m and dips approximately 20 • to the west.

A F I E L D DATA E X A M P L E
To mitigate the effects of displacement currents, Linde & Pedersen (2004) restrict the data set used in quasi-static 2-D inversions with the REBOCC scheme (Siripunvaraporn & Egbert 2000) to frequencies up to 56 kHz. Linde & Pedersen (2004) perform inversions for the TE-mode, TM-mode, TE-and TM-modes together and the determinant of the impedance tensor. By computing synthetic TE-mode, TM-mode and determinant data for a 3-D model and comparing the corresponding 2-D inversions, Pedersen & Engels (2005) show that the inversion of determinant data is less prone to introducing artefacts from 3-D structures off the profile to 2-D inverse models. Furthermore, the inverse model of the determinant data, presented by Pedersen & Engels (2005) has a better data fit than their other models. For the inversion of the RMT data fromÄvrö, this leads us to concentrate on the inversion of determinant data, as the data at both ends of the profile show a high degree of three-dimensionality (Linde & Pedersen 2004). At a few stations, the determinant data of the highest frequencies (160 and 226 kHz) have very small negative phases, which can be indicative of displacement currents (Song et al. 2002). As the rather irregular behaviour of the apparent resistivities at the same stations and frequencies hints at problems with measurement accuracy, we excluded such data points from the inversion.
In the following, we examine the effect of displacement currents, by first considering the inversion of the restricted set of frequencies and then for the full set of frequencies. For each data set, inversions were carried out in both the quasi-static approximation and with displacement  . 8). The resistivity values of borehole KAV01 are taken from the normal-resistivity log presented in Gentzschein et al. (1987). In contrast to models QL and QF, models DL and DF have a more realistic range of resistivities if compared in terms of the range observed in the normal-resistivity log. Furthermore, the resistivity-depth section of model DF at borehole KAV01 is in good agreement with the normal-resistivity log down to a depth of 230 m, and the positions of the seismic reflectors are in good agreement with resistivity contrasts in model DF.
currents. We assumed the relative dielectric permittivity to be r = 6, which, for granite, is in the range between 5 and 8, given by Reynolds (1997). Variation of the permittivity in this range leads to only small differences of the resistivity models for theÄvrö data (not shown).
Our quasi-static determinant model for the lowest frequencies up to 56 kHz (model QL) in Fig. 14(a) resembles the corresponding model by Linde & Pedersen (2004) (their fig. 9d) strongly. We did not include the shallow sea (less than 10 m deep) to the east ofÄvrö as a priori information, as this turned out to be of negligible importance. The central unfractured granite reaches resistivities up to 500 000 m. The conductor at the western end of the profile is interpreted by Linde & Pedersen (2004) as a 150 m wide wet fracture zone, assumed to be related to seismic reflector C of Juhlin & Palm (1999), although the positions of the conductor and reflector are not in very good agreement. The subhorizontal seismic reflector D does not appear to be related to any structure in the resistivity model. The rms misfit of model QL is 1.56.
The inversion with displacement currents for the low-frequency data set gives a model (model DL in Fig. 14b) with a significantly reduced range of resistivities from 300 to 100 000 m. The conductors at 50 and 850 m along the profile appear at greater depth and the boundary between the central resistor and the western conductors is less steep than in the model QL (Fig. 14a). The rms misfit of model DL is 2.03. The quasi-static inversion for the full set of frequencies (model QF in Fig. 14c) leads to a transition into the top of the central resistor that is sharper than in model QL (Fig. 14a). The resistivity of the central resistor is as high as 5 × 10 6 m. The positions of the conductors at profile metres 50 and 850 is very similar to the positions in model QL. The rms misfit of model QF is 3.16.
In comparison to the quasi-static inversions, the inversion with displacement currents for all frequencies gives a model [model DF in Fig. 14(d) with an rms error of 2.60] that shows a less extreme range of resistivities, both at depth and close to the surface. The depth to the conductors at both ends of the profile is about 50 m larger than in the quasi-static models QL and QF in Figs 14(a) and (c), respectively. The model is also in better agreement with the positions of the seismic reflectors. The position of seismic reflector C conforms to an expected boundary between an unfractured resistive granite body and water saturated fractured bedrock. Therefore, we would expect reflector C to represent a boundary of rock units, with different grades of fracturing, rather than a 150 m wide fracture zone, as proposed by Linde & Pedersen (2004). Similarly, reflector D appears to coincide with a subhorizontal boundary of rock units. Furthermore, the model in Fig. 14(d) is in very good agreement with the resistivities of the normal-resistivity log of borehole KAV01 (Gentzschein et al. 1987), down to a depth of 230 m. At greater depth, the model might be more influenced by the smoothness constraint imposed during the inversion than the data. Compared with model DF, model DL (Fig. 14b) deviates from the normal-resistivity log at shallow depth down to 100 m and the positions of Figure 15. Models of the inversion of (a) TE-mode data (model TEDF) and (b) TM-mode data (model TMDF) fromÄvrö for the full set of frequencies allowing for displacement currents. The lines marked by C and D indicate seismic reflectors from Juhlin & Palm (1999). The resistivity values of borehole KAV01 are taken from the normal-resistivity log presented in Gentzschein et al. (1987). The resistivities of the central resistor in model TEDF are as high as 300 000 m. Compared with model DF (Fig. 14d), neither model TEDF nor model TMDF shows similarly good agreement with the positions of the seismic reflectors C and D or the normal-resistivity log.
the seismic reflectors are not as representative as bounds of different rock units. Hence, it appears that the inclusion of high-frequency data is of great importance during the modelling process.
As a verification that the 2-D inverse models of determinant data are less biased by 3-D structures off the profile, the inverse models of TE-mode data (model TEDF) and TM-mode data (model TMDF) are shown in Figs 15(a) and (b), respectively. In both inversions, the full set of frequencies was used and displacement currents were accounted for. The rms fits of 4.56 for the TE-mode model (reached after nine iterations) and 3.67 for the TM-mode model (reached after five iterations) are both significantly higher than that of model DF. The worst data fits of models TEDF and TMDF (not shown) are obtained at the western end of the profile, where strong 3-D effects in the VMT are observed by Linde & Pedersen (2004). In model TEDF, resistivities of the central resistor are as high as 300 000 m. Compared with model DF (Fig. 14d), neither model TEDF nor model TMDF shows similarly good agreement with the positions of the seismic reflectors C and D of Juhlin & Palm (1999) or the normal-resistivity log by Gentzschein et al. (1987).

D I S C U S S I O N A N D C O N C L U S I O N S
We demonstrated the effect of displacement currents on 2-D TM-mode, TE-mode and VMT data, measured with the RMT method at frequencies between 10 and 300 kHz. Forward modelling of subsurfaces with resistivities larger than 1000 m confirms that responses computed in the quasi-static approximation, that is, when displacement currents are neglected, become increasingly inaccurate, with rising frequency. For a homogeneous half-space, both apparent resistivity and phase, computed with displacement currents, decrease from their constant values in the quasi-static approximation, with increasing frequency. At high frequencies, the dielectric effect leads to the occurrence of distinct sign reversals in the real part of the VMT, which are not observed in the quasi-static approximation and might lead to artefactual 2-D or 3-D structures in an interpretation, based on the quasi-static approximation.
The interpretation of high-frequency RMT data with an inverse scheme that operates in the quasi-static approximation will inevitably lead to an inverse model with artefactual structures. As can be seen from the quasi-static interpretation of our synthetic data example in Fig. 11, the resistivities found in this inverse model vary over a larger range than those of the true model (Fig. 9). Typical artefactual structures include conductive near-surface layers, regions of excessively high resistivities next to conductors, as well as conductors that deviate strongly from their true shapes and positions. As only the resistivity distribution is inverted in the scheme presented here, a value for the dielectric permittivity must be chosen before the inversion. The relative dielectric permittivity of bedrock is typically in the range of 5 to 9 (e.g. Reynolds 1997, table 12.3), and a variation in this range does not lead to any important differences in the obtained resistivity models.
Typically, the primary EM field from remote radio transmitters has an angle of incidence that is close to grazing incidence at the measurement site. The assumption of vertically incident plane waves in the modelling code is a limitation, which is of minor importance in many practical situations. Often, a conductive surface layer consisting of, for instance, weathered bedrock or glacial till is present in the area of interest and refracts the incident field towards the vertical due to its relatively low resistivity.
For theÄvrö field data, the inversion that allows for displacement currents and includes high-frequency data produces a model that is in very good agreement with the results of other geophysical methods. The seismic reflectors C and D by Juhlin & Palm (1999) coincide with the boundaries between structures of different conductivity (Fig. 14d). The resistivity depth section of the model at borehole KAV01 matches the normal-resistivity log by Gentzschein et al. (1987) very well, down to a depth of 230 m below which the model might be strongly influenced by the smoothness constraint applied during the inversion. The inverse models computed in the quasi-static approximation (Figs 14a and c) contain artefactual structures, with unrealistically large resistivities, even if only the low-frequency data set is inverted (Fig. 14a).

A C K N O W L E D G M E N T S
The original 2-D MT code by WS at Mahidol University is supported by a research grant from the Thailand Research Fund (TRF: RMU5080025). Nazli Ismail of Uppsala University computed the 3-D forward responses with the X3D code. Dmitry Avdeev, IZMIRAN, Russian Academy of Sciences, spent a large amount of money on the telephone, from Japan, to correct our X3D input files. A discussion with John Weaver, University of Victoria, clarified a few theoretical aspects. We would also like to thank Ute Weckmann, GFZ Potsdam and an anonymous reviewer for their comments, which helped to improve the manuscript.

A P P E N D I X A : F I N I T E -D I F F E R E N C E A P P RO X I M AT I O N
The derivation of the finite-difference approximation (FDA) with the integration method following Hohmann (1987) and Aprea et al. (1997) is particularly instructive. In a finite-difference mesh, a node (i, j) is a corner point of four finite-difference cells, with admittivitieŝ y i−1/2, j−1/2 ,ŷ i+1/2, j−1/2 ,ŷ i−1/2, j+1/2 , andŷ i+1/2, j+1/2 (Fig. A1). The cells have widths y i−1/2 and y i+1/2 and heights z j−1/2 and z j+1/2 . The rectangle A has the centres of the cells as its corner points (Fig. A1). It is assumed that the horizontal electric field component of the TE-mode and the horizontal magnetic field component of the TM-mode at node (i, j) are E i, j x and H i, j x , respectively, and that the magnetic permeability is equal to its vacuum value μ 0 . Nodes along the boundary of the finite-difference mesh are called boundary nodes. All other nodes are called inner nodes. Finite-difference equations for the TM-and TE-mode are obtained by integrating eqs 14 and 13, respectively, over the surface of A. This surface integral is then transformed to a contour integral around the perimeter ∂ A of the surface A with Gauss' Theorem.
Hence, for the TE-mode, where n is an outward unit normal vector on the edges of A. The part of A, for instance, which is entirely situated to the upper left-hand side of node (i, j), contributes to the surface integral in eq. (A1) with For the upper edge of the rectangle A, for instance, the integral around the perimeter ∂ A can be approximated as where ( y i+1/2 + y i−1/2 )/2 is the length of the perimeter ∂ A on the upper edge of A and n, which equals (0, 0, −1) on the upper edge of A, collects the vertical component of ∇ E x , that is, )/ z j−1/2 , multiplied by −1. In total, this leads to the following approximations (cf . Aprea et al. 1997): Figure A1. Finite-difference mesh around node (i, j). The four surrounding cells have admittivitiesŷ i−1/2, j−1/2 ,ŷ i+1/2, j−1/2 ,ŷ i−1/2, j+1/2 andŷ i+1/2, j+1/2 . The heights and widths of the surrounding cells are z j−1/2 and z j+1/2 and y i−1/2 and y i+1/2 , respectively. The rectangle A has its corner points at the centres of the cells.
After obtaining j i, j y from the above equations, the two one-sided values E i−, j y and E i+, j y of the electric field component E y can be computed. In contrast to eq. (24) in Weaver et al. (1985), the current density j i, j y must be divided by the left-and right-hand sided vertically averaged admittivities (eqs A22 and A21, respectively) to obtain E i−, j y and E i+, j y , respectively, that is, The use of averaged admittivities is motivated by considering the integrated current I y = x−z− plane j y dxdz through any surface y = const.
To assign a unique value to E i, j y , the current density j i, j y is typically divided by the average admittivityŷ avg i, j given in eq. (A7), that is, y avg i, j E i, j y = j i, j y .
where the computation of u k as a solution of the pseudo-forward problem K T u k =c k exploits the reciprocity of the forward problem (i.e. the symmetry of the forward system matrix K and its inverse). The computation of u k for k = 1, . . . , N is significantly faster than the computation of K −1 (− ∂K ∂m l x + ∂s ∂m l ) for l = 1, . . ., M in the case of typical 2-D problems, where the number of model parameters exceeds the number of data, that is, where M > N .
A model parameter m l = ρ i+1/2, j+1/2 is connected to the cell to the lower right-hand side of a node (i, j) through l = j − N za + (i − 1) * N zb . Similarly, each data index k is connected to a single surface node (i s , j s = N za + 1) for a given frequency, where i s and j s indicate the node at which a certain receiver station is located. were used. The quantityŷ avg i, j is given according to eq. (A7). The derivative of the system matrix of the forward problem K w.r.t. a single model parameter m l in eq. (C7) results in a matrix ∂K/∂m l that has only four rows with non-zero entries. The parameter m l = ρ i+1/2, j+1/2 enters into the rows of K that correspond to the central nodes (i, j), (i, j + 1), (i + 1, j) and (i + 1, j + 1) (cf. Fig. A1), that is, into rows number The computation of ∂K/∂m l is further simplified by the symmetry of K. For a central node (i, j) the coefficient of the EM field component at its right-hand side node is the same as the coefficient of the EM field component at the left-hand side node of its neighbouring central node (i + 1, j). Similarly, for a central node (i, j) the coefficient of the EM field component at its lower node is the same as the coefficient of the EM field component at the upper node of its neighbouring central node (i, j + 1). Furthermore, as the left-and right-hand coefficients contain vertically averaged inverse admittivities and the lower and upper coefficients contain horizontally averaged inverse admittivities, the derivative of the coefficient of the right-hand node of the central node (i, j) w.r.t. ρ i+1/2, j+1/2 equals the derivative of the coefficient of the right-hand node of its neighbouring central node (i, j + 1) w.r.t. ρ i+1/2, j+1/2 . Similar rules are valid for the coefficients of left-hand, upper and lower nodes at correspondingly neighbouring nodes. Hence, ∂K(iul, iul + (N z − 1)) ∂ρ i+1/2, j+1/2 = ∂K(idl, idl + (N z − 1)) ∂ρ i+1/2, j+1/2 = ∂K(iur, iur − (N z − 1)) ∂ρ i+1/2, j+1/2 = ∂K(idr, idr − (N z − 1)) ∂ρ i+1/2, j+1/2 = 2 z j+1/2 y i+1/2

C2 TE-mode
In the TE-mode, H y is expressed through b k and E x according to eq. (A17) and a k is zero except for the kth entry, which is 1. Hence, The Authors, GJI, 175, 486-514 Journal compilation C 2008 RAS