# Jointly reconstructing ground motion and resistivity for ERT-based slope stability monitoring

Jointly reconstructing ground motion and resistivity for ERT-based slope stability monitoring Summary Electrical resistivity tomography (ERT) is increasingly being used to investigate unstable slopes and monitor the hydrogeological processes within. But movement of electrodes or incorrect placement of electrodes with respect to an assumed model can introduce significant resistivity artefacts into the reconstruction. In this work, we demonstrate a joint resistivity and electrode movement reconstruction algorithm within an iterative Gauss–Newton framework. We apply this to ERT monitoring data from an active slow-moving landslide in the UK. Results show fewer resistivity artefacts and suggest that electrode movement and resistivity can be reconstructed at the same time under certain conditions. A new 2.5-D formulation for the electrode position Jacobian is developed and is shown to give accurate numerical solutions when compared to the adjoint method on 3-D models. On large finite element meshes, the calculation time of the newly developed approach was also proven to be orders of magnitude faster than the 3-D adjoint method and addressed modelling errors in the 2-D perturbation and adjoint electrode position Jacobian. Hydrogeophysics, Electrical resistivity tomography (ERT), Inverse theory, Tomography 1 INTRODUCTION Electrical resistivity tomography (ERT) is a geophysical method where current is injected between pairs of electrodes attached to the surface of the earth while the voltage differences are measured between other electrode pairs. From this data, one may reconstruct a resistivity distribution for the near surface earth that best fits the available data. For ERT, relatively small electrode movements or boundary modelling errors can lead to reconstructions with resistivity artefacts so severe that the resulting image is difficult to interpret (Zhou & Dhalin 2003). Similar reconstruction artefacts have been observed for biomedical electrical impedance tomography (EIT) where the fundamental mathematical problem is identical to ERT (Adler et al. 1996, 2011; Boyle & Adler 2011). Intermittent electrode movement is expected during long-term monitoring of an active landslide site. One could re-survey the electrode locations after each movement, but this would be time consuming and expensive, particularly for remote locations. Changes in landslide topography are frequently accompanied by changes in the resistivity distribution of the slope: changes in soil properties, such as water saturation, affect slope stability and resistivity. A resistivity reconstruction that does not account for electrode movement, when movements have occurred, will tend to fit the data poorly. Attempts to force fit a resistivity distribution will, in most cases, lead to images with severe resistivity artefacts which, without further information, could be misinterpreted. Simultaneously adjusting electrode position and resistivity can allow for better data fit, reduced resistivity artefacts, and indications of electrode movement. The methods described here separate the confounding factors, electrode movement and resistivity changes, which have long hindered ERT landslide monitoring attempts. Previous work on electrode movement has largely focused on 2-D electrode position changes in the plane of the electrodes. In 2-D, changes in boundary shape that are conformal do not result in artefacts (Boyle et al. 2012a). Modelling in 3-D is important because it reflects the finite size of the electrodes and the resulting current flow. For 3-D reconstructions, conformal deformations are limited to rotation, scaling and translation which would normally be identified through an appropriate site survey. 3-D resistivity reconstructions can be prohibitive to calculate for a given level of fidelity and numerical convergence. The so-called 2.5-D method (Dey & Morrison 1979) combines a 2-D parametrization of the resistivity and electrode positions with a 3-D model of current flow in the ground. The method can enable significant reductions in the computational requirements for meshing and current density calculations while producing accurate results with respect to equivalent 3-D models. Electrode movement and resistivity have been reconstructed for biomedical problems where changes both electrode movement (±1 per cent of electrode spacing) and resistivity (±50 per cent) have been relatively mild, enabling Gauss–Newton difference solutions (Soleimani et al. 2006; Jehl et al. 2015), and for large deformations on cylindrical domains with surface-normal deformations (Dardé et al. 2013a,b; Hyvönen et al. 2014) in simulation and on tank data. For geophysics, the resistivity contrasts are generally strong, varying by orders of magnitude. These strong contrasts combined with large electrode movements, much greater than 1 per cent of electrode spacing, necessitate an iterative solution because the combined effects of the electrode movement give highly nonlinear interactions. Similar approaches are also currently being developed by other researchers (Kim et al. 2014; Wagner et al. 2015; Loke et al. 2015, 2016, 2017). Results on field data have been somewhat mixed and data dependent. One approach to account for electrode position modelling errors is to alternate between an electrode position update and resistivity update (Wilkinson et al. 2010, 2016). Such an orthogonal greedy optimization approach may be slow to converge, or even fail to converge, because descent over the optimization surface is restricted to movement parallel to the axes (Cormen et al. 1990). Ideally, both electrode movement and resistivity should be reconstructed simultaneously so that updates can descend in the globally optimal direction across the objective function. In this work, joint electrode position and resistivity reconstruction methods were applied, in 2.5-D, to two data sets (Figs 1 and 2, line#1 and line #5) from an active landslide in the UK. An efficient and accurate method for calculating the electrode position Jacobian on a 2.5-D resistivity model was developed. This work is motivated by developments in Boyle & Adler (2010); Boyle (2010) for 2-D electrode movement, and Boyle et al. (2014); Boyle (2016); Boyle & Adler (2016) where this data set and preliminary results were presented. A 2-D cross-sectional Finite Element Method (FEM) model of the local slope topology was built, with electrodes modelled using the Complete Electrode Model (CEM) (Somersalo et al. 1992; Rücker & Günther 2011). We reconstruct electrode movement and resistivity in an iterative Gauss–Newton algorithm, showing that under certain conditions both can be simultaneously reconstructed. Figure 1. View largeDownload slide Slope failures at Hollin Hill; (a) rotational failures near the top of the slope above line#5, June 2015, (b) earth flow at the toe of a landslide where line#5 runs through mid-slope with electrodes throughout the slipped material, June 2015, (c) satellite image of the hillside (2016), showing four landslide ‘lobes’, five lines of 32 electrodes as of February 2014, and base station location.[(c) Satellite imagery ©2016 DigitalGlobe, Getmapping plc, Infoterra Ltd Bluesky.] Figure 1. View largeDownload slide Slope failures at Hollin Hill; (a) rotational failures near the top of the slope above line#5, June 2015, (b) earth flow at the toe of a landslide where line#5 runs through mid-slope with electrodes throughout the slipped material, June 2015, (c) satellite image of the hillside (2016), showing four landslide ‘lobes’, five lines of 32 electrodes as of February 2014, and base station location.[(c) Satellite imagery ©2016 DigitalGlobe, Getmapping plc, Infoterra Ltd Bluesky.] Figure 2. View largeDownload slide Electrode locations; (a) electrode locations for five lines of thirty-two electrodes each, line#1 (blue) to line#5 (green) as of February 2014, (b) where electrodes were 10 mm × 170 mm spikes of stainless steel selected for its conductivity, low cost and corrosion resistance. Figure 2. View largeDownload slide Electrode locations; (a) electrode locations for five lines of thirty-two electrodes each, line#1 (blue) to line#5 (green) as of February 2014, (b) where electrodes were 10 mm × 170 mm spikes of stainless steel selected for its conductivity, low cost and corrosion resistance. 2 BACKGROUND Resistivity imaging has been used in geophysical investigations of the behaviour and precursors of landslides and failure surfaces (Perrone et al. 2004; Lapenna et al. 2005; Lebourg et al. 2005; Jongmans & Garambois 2007; Naudetb et al. 2008; Sass et al. 2008). The technique is attractive because resistivity is strongly dependent on water saturation, fracturing, clay content and weathering which are all key factors in slope stability (Piegaria et al. 2009). Slopes may be monitored over time to observe changes in these key parameters using automated systems to collect and analyse data on a daily basis (Kuras et al. 2009; Lebourg et al. 2010; Supper et al. 2014). Difference images may show immediate changes in water saturation (Suzuki & Higashi 2001; Friedel et al. 2006; Jomard et al. 2007) but are limited in their ability to perform long-term monitoring due to background resistivity changes and electrode movements. As illustrated in Fig. 3, small amounts of electrode movement may introduce significant artefacts (Zhou & Dhalin 2003; Oldenborger et al. 2005; Wilkinson et al. 2008). These artefacts may be reduced by accounting for changes in electrode position (Uhlemann et al. 2017). Figure 3. View largeDownload slide Electrode movement artefacts; simulated reconstructions with a conductive and insulating target (rectangular and square outlines, respectively), each with two electrodes (electrode #2 and #12 of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (a) 0 per cent, (b) 5 per cent and (c) 25 per cent of electrode spacing on a 2-D half-space reconstruction (40 dB Signal-to-Noise Ratio (SNR), λ = 0.01, Laplace regularization, Wenner measurement pattern). Single or well separated electrode location errors introduce characteristic ‘ringing’ artefacts that can overwhelm resistivity-based information. [Reproduced, in part, from Boyle et al. (2017), fig. 2., with author permission.] Figure 3. View largeDownload slide Electrode movement artefacts; simulated reconstructions with a conductive and insulating target (rectangular and square outlines, respectively), each with two electrodes (electrode #2 and #12 of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (a) 0 per cent, (b) 5 per cent and (c) 25 per cent of electrode spacing on a 2-D half-space reconstruction (40 dB Signal-to-Noise Ratio (SNR), λ = 0.01, Laplace regularization, Wenner measurement pattern). Single or well separated electrode location errors introduce characteristic ‘ringing’ artefacts that can overwhelm resistivity-based information. [Reproduced, in part, from Boyle et al. (2017), fig. 2., with author permission.] An active landslide was identified in North Yorkshire, UK (54°06΄39.2″N, 0°57΄34.9″W) and has been monitored since 2008 (Chambers et al. 2011; Wilkinson et al. 2016). The landslide is a very slow to slow moving composite multiple earth slide-earth flow (Uhlemann et al. 2016b), where a central portion of the slope has moved downhill by up to 3.5 m a year in some instances (Fig. 1). The slope itself exposes four formations: the Dogger Formation (DF), Whitby Mudstone Formation (WMF), Staithes Sandstone Formation (SSF) and Redcar Mudstone Formation (RMF), from top to bottom. The interfaces between sedimentary layers lie horizontally, with a gentle 5° dip to the North, determined through comparison of material interfaces at surrounding exposed slopes in the region (Chambers et al. 2011). The WMF, as the name implies, is a mudstone clay-based rock that is highly weathered and prone to movement during peak water saturation periods at Hollin Hill, occurring annually in the winter through early spring wet-season. The underlying SSF and unweathered RMF are structurally more competent, and landsliding is postulated to occur within the WMF (Uhlemann et al. 2016b). The slope lies at an average angle of 14° over a change in elevation of 40 m. The landslide movement is determined by translational movements upslope of the WMF–SSF interface. These lead to a loss of support of the local slope of the back-scarp, causing rotational failure. The mass accumulation at the WMF–SSF interface, driven by these translational movements, and elevated pore pressures then cause flow movements which form lobes (Uhlemann et al. 2016a,b). A grid of five rows of thirty-two permanently installed electrodes travelled along with these movements, shifting their positions relative to each other (Fig. 2). An automated impedance imaging survey was executed bi-daily and data were transmitted to a remote office for storage and analysis. In 2008–2009, a middle section of line#1 exhibited a translational failure with movements of up to 1.6 m. In 2013–2014, upper and middle segments of line#5 had rotational and translational failures of up to 3.5 m. The dipole-dipole measurement protocol used for line#1 and line#5 are visualized in Figs 4 and 5, showing the sequence of quadrupolar measurements with stimulus dipoles in red and measurement dipoles in blue. A single difference measurement is captured as one row of the protocol in the figure. In the adjacent vertical strip chart, the horizontally aligned measurements at the initial time (green) is contrasted with the homogeneous resistivity estimate (red) of what those measurements would be, and the change in measurements from initial to final measurement (blue). The rightmost strip chart shows the estimated error for each measurement as estimated from differences between reciprocal pairs of data for the initial measurements (green) and final measurements (blue). The figure serves to illustrate three challenges with this data set. First, the measurements do not fit, or nearly fit, a homogeneous resistivity model. Second, the change between the initial and final measurements (resistivity and movement) is as much as that caused by the inhomogeneous resistivity distribution in the initial measurements (resistivity only) for some data. Third, the measurement error varies over time so that the initial and final measurements have different associated error estimates. The sequence of measurements for line#5 (Fig. 5) started at the top of the slope and ran to the base, but is equivalent to that shown for line#1 (Fig. 4). The final measurements for line#5 have significantly worse reciprocal errors, but are the best of many data sets examined for a post-movement data set on line#5. Figure 4. View largeDownload slide Dipole–dipole measurement protocol for line#1; March 2008 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 32.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 4. View largeDownload slide Dipole–dipole measurement protocol for line#1; March 2008 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 32.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 5. View largeDownload slide Dipole–dipole measurement protocol for line#5; February 2013 and February 2014 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 26.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 5. View largeDownload slide Dipole–dipole measurement protocol for line#5; February 2013 and February 2014 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 26.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. 3 PRIOR RECONSTRUCTIONS AT HOLLIN HILL Electrode movement has been previously reconstructed for line#1 data (2008/2009) along with separate resistivity sections (Wilkinson et al. 2010). The algorithm used in that instance achieved estimates of within 0.2 m (4 per cent of electrode spacing) of the electrode’s true positions as measured by Real-Time Kinematic Global Positioning System (RTK GPS).1 An initial resistivity-only reconstruction with known electrode positions gave a plausible distribution that was in good agreement with available geological evidence: borehole and auger data, evaluation of local geology, aerial LiDAR, differential GPS, lab correlation of representative samples to measured conductivities, piezometric pore pressure measurements, and in-situ rainfall and temperature records (Fig. 6; Chambers et al. 2011; Merritt et al. 2014). Resistivity reconstructions of data collected after movement exhibited artefacts. These artefacts were reduced when reconstructed movements were incorporated. The electrode movement was penalized in the upslope direction. The position Jacobian was estimated based on an analytic half-space model with a homogeneous resistivity assigned to each group of measurements based on electrode separation. The electrode movement was then reconstructed by minimizing   \begin{eqnarray} \arg \min \sqrt{\sum _i |e_i|^2} + \alpha \sum _j |m_j| + \beta \sum _j \theta (m_j)|m_j| \end{eqnarray} (1)for weighting factors α = 0.06 m−1, β = 0.32 m−1, Heaviside step function θ, movement mj, and misfit ei = rb − ra as the difference between predicted rb and observed ra apparent resistivity ratios. The apparent resistivity ratio was calculated as the ratio of the analytic half-space models (9) before and after movement   \begin{eqnarray} r = \frac{\rho ^{\prime }}{\rho } \left( \frac{ \frac{1}{AM^{\prime }} - \frac{1}{BM^{\prime }} - \frac{1}{AN^{\prime }} - \frac{1}{BN^{\prime }}}{\frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} - \frac{1}{BN}} \right) \end{eqnarray} (2)for homogeneous resistivity ρ, electrodes spaced AM, BM, AN, BN, and the updated locations and resistivity after movement indicated by primed ΄ variables. Dipole-dipole data for measurements n = 1 were discarded, as they were judged to be too dependent on transverse movements which were not reconstructed. A similar approach was applied (Wilkinson et al. 2016), to reconstruct 2-D surface xy-movements for the whole electrode grid by allowing for transverse movements through an additional weighting term   \begin{eqnarray} \arg \min &&\sum _i e_i^2 + \alpha \sum _j ||m_j|| + \beta \sum _j \theta (m_j^{(y)})|m_j^{(y)}| \nonumber\\ &&\quad+ \, \gamma \sum _j \theta (m_j^{(x)})|m_j^{(x)}| \end{eqnarray} (3)facilitating balancing the sensitivities of transverse $$m_j^{(y)}$$ and longitudinal $$m_j^{(x)}$$ movements by adjusting the weighting ratio β/γ. Figure 6. View largeDownload slide [From Wilkinson et al. (2010), figure 2, x, z-coordinates corrected:] 2-D resistivity image inverted from the baseline data set (2008 March). The inferred boundaries between the Whitby (WMF), Staithes (SSF) and Redcar (RMF) formations are shown by dotted black lines. Stratigraphic logs of boreholes are shown in grey scale. The main scarp and slipped WMF material are indicated by the black arrows. [Reproduced from Wilkinson et al. (2010), fig. 2, for comparison with Fig. 11(a) in this work. The RES2DINV reconstruction region, elsewhere in this paper referred to as the ‘RES2DINV outline,’ is selected by the software based on a heuristic pseudo-section method (Loke 2017).] Figure 6. View largeDownload slide [From Wilkinson et al. (2010), figure 2, x, z-coordinates corrected:] 2-D resistivity image inverted from the baseline data set (2008 March). The inferred boundaries between the Whitby (WMF), Staithes (SSF) and Redcar (RMF) formations are shown by dotted black lines. Stratigraphic logs of boreholes are shown in grey scale. The main scarp and slipped WMF material are indicated by the black arrows. [Reproduced from Wilkinson et al. (2010), fig. 2, for comparison with Fig. 11(a) in this work. The RES2DINV reconstruction region, elsewhere in this paper referred to as the ‘RES2DINV outline,’ is selected by the software based on a heuristic pseudo-section method (Loke 2017).] 4 METHOD While the work of Wilkinson et al. (2010, 2016), described in the previous section, focused on a sequential inversion procedure, here we simultaneously reconstruct resistivity and electrode position in a Gauss–Newton iterative framework. A 2-D model of the slope profile was constructed with independent parameters for resistivity and electrode position (Fig. 7). The ℓ2-norm of the data discrepancy and regularization terms was minimized by balancing the sensitivity of the resistivity parameters ρ and electrode movement x against the regularization terms   \begin{eqnarray} (\hat{\mathbf {\rho }},\hat{{\mathbf {x}}}) &=& \arg \min _{\mathbf {\rho },{\mathbf {x}}} \ ||\mathcal {F}(\mathbf {\rho },{\mathbf {x}}) - {\mathbf {m}}||_{\mathbf {W}} \nonumber \\ &&\quad + \, ||\lambda \beta {\mathbf {R}}_\mathbf {\rho }\log _{10}(\mathbf {\rho }_\star /\mathbf {\rho })||_{2} \nonumber \\ && \quad + \, || \lambda \eta {\mathbf {R}}_{{\mathbf {x}}} ({\mathbf {x}}- {\mathbf {x}}_\star )||_{2}. \end{eqnarray} (4)The optimal solution $$(\hat{\mathbf {\rho }},\hat{{\mathbf {x}}})$$ minimizes the data discrepancy between a forward model $$\mathcal {F}$$ and measurements m combined with some regularization Rρ and Rx acting to smooth an otherwise ill-posed and ill-conditioned inverse problem. The forward model was the parametrized 2.5-D model (Fig. 7) of resistivity ρ and longitudinal2 electrode position x and producing an estimate of expected measurements given the measurement protocol. Measurement error was weighted W based on estimates of measurement reliability. The regularization terms penalized changes in resistivity and electrode position from prior estimates (ρ⋆, x⋆). The relative sensitivity of the two types of parameters, resistivity and electrode movement, were balanced by adjusting the ratio of the scalar β/η. The overall strength of the regularization was adjusted by scaling both terms by λ. The resistivity was solved under a positive conductivity σ constraint by converting to inverse log units log 10σ = log 101/ρ. Figure 7. View largeDownload slide The 2.5-D forward model FEM parametrization; right-to-left expanding from the region surrounding a single electrode, to the scale of the electrode array, to the scale of the region surrounding the electrode array, the forward model is parametrized for electrode position x and resistivity ρ, resistive regions selected to demonstrate mesh structure. Figure 7. View largeDownload slide The 2.5-D forward model FEM parametrization; right-to-left expanding from the region surrounding a single electrode, to the scale of the electrode array, to the scale of the region surrounding the electrode array, the forward model is parametrized for electrode position x and resistivity ρ, resistive regions selected to demonstrate mesh structure. The objective function (4) was solved using the well-known iterative Gauss–Newton approach (Nocedal & Wright 1999,  section 10.2). The Gauss–Newton approach starts from an initial estimate (ρ0, x0), estimates the local slope of the objective function as the Jacobian (Jσ, Jx) to determine a new search direction (δρ, δx), and then performs an approximate line search in that direction to estimate an optimal step length α. At each iteration, the parameters were updated   \begin{eqnarray} (\delta \mathbf {\rho }, \delta {\mathbf {x}}) = -({{\bf J}}^{{\sf T}} {\mathbf {W}}{{\bf J}}+ \lambda ^2 {\bf Q})^{-1}({{\bf J}}^{{\sf T}} {\mathbf {W}}{\mathbf {b}}\, + \lambda ^2 {\bf Q} {\bf c}) \end{eqnarray} (5)  \begin{eqnarray} &&{(\mathbf {\rho }_{n+1}, {\mathbf {x}}_{n+1}) = (\mathbf {\rho }_n, {\mathbf {x}}_n) + \alpha (\delta \mathbf {\rho }, \delta {\mathbf {x}})} \\ &&{{\rm for }\, \,{\mathbf {b}}= \mathcal {F}(\mathbf {\rho }_n,{\mathbf {x}}_n) - {\mathbf {m}}}\nonumber \\ &&{{\bf c} = \left[ \begin{array}{c}\log _{10} 1/\mathbf {\rho }_n - \log _{10} 1/\mathbf {\rho }_\star \\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array} \right] = \left[ \begin{array}{c}\log _{10} \mathbf {\rho }_\star /\mathbf {\rho }_n\nonumber\\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array}\right]}\nonumber \\ &&{{{\bf J}}= \left[ \begin{array}{c}\frac{-1}{\mathbf {\rho }_n \ln (10)} \, {{\bf J}}_{\sigma _n} \\ {{\bf J}}_{\mathbf {x}}\end{array} \right]} \nonumber \\ &&{{\mathbf {R}}= \left[ \begin{array}{c{@}{\quad}c} \beta {\mathbf {R}}_\mathbf {\rho } & 0 \\ 0 & \eta {\mathbf {R}}_{\mathbf {x}}\end{array}\right] \qquad {\rm and \, \, \, {\mathbf {R}}^{{\sf T}} {\mathbf {R}}= {\bf Q}}}\nonumber \end{eqnarray} (6)based on the data discrepancy b and the distance from the prior estimate $$\bf c$$ in combination with the regularization R, the slope of the objective function J and measurement weighting W. This formulation agrees with that of Boyle et al. (2017), but is modified to address a log resistivity parametrization. In contrast to a typical resistivity-only inversion, the reconstruction parameters, Jacobian and regularization R have been extended to incorporate the new electrode position parameters. The resistivity Jacobian has been calculated for conductivity $${{\bf J}}_{\sigma _n}$$ and then scaled, which is exactly equivalent to calculating the Jacobian on the log of resistivity. The movement Jacobian was found to be sensitive to resistivity changes in the small elements adjacent to the electrodes. To address this sensitivity, the log conductivity regularization combined a smoothing prior near electrodes with Tikhonov regularization away from the electrodes   \begin{eqnarray} {\mathbf {R}}_\rho = {\mathbf {I}}+ \nu {\mathbf {R}}_{\mathfrak {L}} \end{eqnarray} (7)  \begin{eqnarray} \nu = \exp -\, \frac{|x_e - x_\ell |}{|\bar{x}_{\ell \ell }|} \end{eqnarray} (8)where the Laplace smoothing $${\mathbf {R}}_{\mathfrak {L}}$$ was scaled ν by the distance between each FEM element centre xe and the nearest electrode xℓ, scaled by the average distance between electrodes $$\bar{x}_{\ell \ell }$$. This prior encourages small changes from the expected resistivity in regions with little information. In regions near the electrodes, changes will be pushed towards the spaces between electrodes rather than directly under the electrodes, as well as encouraging smooth transitions in resistivity near the electrodes. The regularization for movement was the Tikhonov prior Rx = I. In principle, there are correlated changes between resistivity and electrode movement (Kim et al. 2014) which may be partially accounted for by setting the off-diagonal blocks of the regularization matrix to non-zero values, but in practice these were not characterized and in the absence of a better guess were set to zero. The forward model was constructed as a 2-D cross-section based on the original electrode locations and the mesh was then perturbed by PCHIP3 interpolation (Carlson & Fritsch 1989) for electrode displacements. Forward modelled measurements and the Jacobians were calculated using the 2.5-D method (Dey & Morrison 1979). We have made use of the log conductivity to restrict the resistivity reconstruction to physically meaningful positive values. Experiments with the log movement constraint, to restrict electrodes to downslope movement, resulted in similar reconstructions to the ones presented here which used unrestricted electrode movement. For the log movement parametrization, the behaviour at each iteration was different to an unscaled movement due to the structure of the movements in this data set. Because each electrode that moved had a different magnitude of movement, the log movement reconstruction tended to solve for each electrode’s reconstructed displacement separately: one electrode per iteration. The apparent single electrode updates were actually an artefact of the log scaling, where smaller movements were reduced by orders of magnitude, so as to be inconsequential. Once the largest electrode placement error had been corrected, the next largest error would be addressed. We feel that this highlights the importance of careful selection of the reconstruction parametrization. It is possible that some Fourier decomposition or other basis of electrode movement with appropriate regularization might achieve greater reconstruction accuracy without artificially fixing any single electrode’s location or degrees of freedom for movement. 5 2.5-D POSITION JACOBIAN The 2-D electrode position Jacobian suffers from significant errors (Fig. 8), when compared to data from a 3-D model, which makes it inappropriate for 3-D problems. A 3-D electrode position Jacobian becomes prohibitively expensive to calculate as the mesh density grows. A 2.5-D approach offers a compromise by restricting sensitivity parametrization and electrode positions to the plane, while achieving high fidelity to equivalent 3-D models, at a fixed multiple of the 2-D computational effort. Figure 8. View largeDownload slide Jacobian sensitivity $${\rm diag}\,\,({{\bf J}}_x^{{\sf T}} {{\bf J}}_x)^\frac{1}{2}$$ on a 16-electrode, homogeneous (σ = 1) half-space model; 2-D rank-one electrode position Jacobian (Gómez-Laberge & Adler 2008) shows orders-of-magnitude error in sensitivity estimate, while 3-D analytic (9) and rank-one estimates (Gómez-Laberge & Adler 2008) are in close agreement with the 2.5-D (13) estimate. Figure 8. View largeDownload slide Jacobian sensitivity $${\rm diag}\,\,({{\bf J}}_x^{{\sf T}} {{\bf J}}_x)^\frac{1}{2}$$ on a 16-electrode, homogeneous (σ = 1) half-space model; 2-D rank-one electrode position Jacobian (Gómez-Laberge & Adler 2008) shows orders-of-magnitude error in sensitivity estimate, while 3-D analytic (9) and rank-one estimates (Gómez-Laberge & Adler 2008) are in close agreement with the 2.5-D (13) estimate. Two alternate methods were evaluated before developing the 2.5-D position Jacobian: a 2.5-D perturbation method which was relatively slow, and an analytic model of movement which was restricted to a homogeneous resistivity and the Point Electrode Model (PEM). Motivated by the short comings of these two methods, we develop the 2.5-D position Jacobian which is efficient and accounts for resistivity variation in the model using the CEM. The perturbation method used the underlying 2.5-D forward simulations of the full FEM resistivity model using the CEM. These movement perturbation calculations proved to be prohibitively slow because a mesh perturbation resulted in recalculations of the system matrices and a new inversion of that system matrix. The cost grows with number of electrodes and movement dimensions, so that for n = 32 electrodes, estimated in d = 1 dimensions, nd + 1 forward simulations were required at each iteration of the Gauss–Newton inverse solution. A line search typically required three to six forward simulations, meaning that the perturbation Jacobian required far more time to calculate than the rest of each iteration. When FEM mesh density was increased sufficiently to achieve good estimates of electrode position changes, the calculations took an unreasonable amount of time. Low mesh densities sped up the calculations but exhibited significant errors when compared to a 3-D perturbation solution at sufficient mesh density. An alternate solution was implemented by adapting the half-space analytic PEM forward model. For a half-space, the potential difference V measured over a homogeneous resistivity ρ with current I driven on the stimulus electrodes, is given by   \begin{eqnarray} V = \frac{\rho I}{2 \pi } \left( \frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} + \frac{1}{BN} \right) \end{eqnarray} (9)where each distance AM, BM, AN, BN is between a stimulus electrode and a measurement electrode. The model may be applied for any arbitrary pair-wise electrode placement. A position Jacobian may be constructed by applying the electrode movement perturbation. A next logical step would be to take the derivative of (9) with respect to electrode position and build up a solution specific approximate block-wise model of resistivity away from electrodes (Wilkinson et al. 2010), though this was not implemented in this work. Electrode positions were captured from the FEM model. A homogeneous resistivity was assigned based on the average resistivity of the current FEM model. The electrode position Jacobian produced by the half-space analytic perturbation method was compared to the 2.5-D perturbation Jacobian under homogeneous conditions. Using the modified half-space analytic perturbations was much faster to calculate than the 2.5-D perturbation method and reasonably accurate: the effect of topography was found to be somewhat accounted for. The loss of accuracy due to changing electrode models (CEM to PEM) and using a homogeneous resistivity were not so disruptive as to change signs in the position Jacobian although magnitudes were inaccurate. Motivated by the efficiency of the 2-D position Jacobian of Gómez-Laberge & Adler (2008), and the errors introduced by using a 2-D electrode position Jacobian, the 2.5-D position Jacobian was developed as the corollary of the 2.5-D conductivity Jacobian. In general, the 2.5-D forward solver is a well-known technique and is commonly used in geophysics ERT applications. An approximately half-space geometry, and a resistivity that is nearly uniform along one axis, fit well with a 2.5-D model, and occur naturally in many geological settings. By adapting the adjoint, or ‘standard method,’ of calculating the conductivity Jacobian by a rank-one update, to the 2.5-D technique, resistivity may be efficiently reconstructed. To reconstruct electrode movement, we desire a similar 2.5-D implementation for the electrode position Jacobian. In the following, we outline the 2.5-D conductivity Jacobian and present a new derivation for the 2.5-D electrode position Jacobian. We use the formulation of the 2-D conductivity Jacobian from Boyle et al. (2017) as a basis for these developments. The 2-D and 2.5-D conductivity Jacobians Jσ were calculated   \begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} {\mathbf {S}}\frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}{\mathbf {X}} \end{eqnarray} (10)  \begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} ({\mathbf {S}}+ k^2 {\mathbf {T}}) \frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray} (11)for measurement selection $${\mathcal {T}}$$, system matrix A, mesh connectivity matrix C, mesh shape functions S, a conductivity change ∂D/∂σn, and the nodal voltages X = A−1B for stimulus B over an electrode modelled as a shunt in the y-direction when the 2-D FEM is meshed over the x–z plane. The system matrix $${\mathbf {A}}= {\mathbf {C}}^{{\sf T}} {\mathbf {S}}{\mathbf {D}}{\mathbf {C}}$$ is assembled from a connectivity matrix C mapping global node numbers to element-local node numbers, the element shape functions S and the conductivity D per element. For the 2.5-D position Jacobian, the system matrix Ak is specific to the spatial frequency k, as are the nodal voltages $${\mathbf {X}}_k = {\mathbf {A}}_k^{-1} {\mathbf {B}}$$. A perturbation node is selected at row u, column v, affecting a linear interpolatory shape function E. The 2.5-D position Jacobian may be calculated as an extension of the 2-D Jacobian, in an analogous way to (10), (11), as   \begin{eqnarray} {{\bf J}}_{x,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial {\mathbf {S}}}{\partial x_n} {\mathbf {D}}{\mathbf {C}}{\mathbf {X}} \end{eqnarray} (12)  \begin{eqnarray} {{\bf J}}_{x,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial ({\mathbf {S}}+ k^2 {\mathbf {T}})}{\partial x_n} {\mathbf {D}}{\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray} (13)where the 2-D position Jacobian may be efficiently calculated using the rank-one update for the conductivity Jacobian (Gómez-Laberge & Adler 2008) with some new terms. Again, based on the 2-D formulation from Boyle et al. (2017), the element shape functions for element (e) may be summarized as the element shape matrix   \begin{eqnarray} {\mathbf {S}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|}{\mathbf {E}}_{\setminus 1}^{{\sf T}} {\mathbf {E}}_{\setminus 1} \end{eqnarray} (14)for a shape matrix E and a row-reduced version E∖1 where the top row of the matrix is removed. The shape matrix is distorted by having its nodes perturbed leading to the first-order estimate of element deformation   \begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} &=& \frac{1}{2}\, \Bigg ( \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} {\mathbf {E}}_{\setminus 1}^{{{\sf T}} } {\mathbf {E}}_{\setminus 1} \, \nonumber\\ &&\qquad + \frac{1}{|\det {\mathbf {E}}|} \bigg (\frac{\partial {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }}{\partial x_n} {\mathbf {E}}_{\setminus 1} + {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }\frac{\partial {\mathbf {E}}_{\setminus 1}}{\partial x_n} \bigg ) \Bigg ) \end{eqnarray} (15)where xn refers to a global node numbered n and affects all elements e connected to that node. The local shape functions of each element, for first-order interpolatory shape functions on a 2-D mesh, are   (16)for a triangular element (blue) with three nodes p1, p2, p3. To calculate the partial derivatives of the first-order interpolatory shape functions, we make use of the matrix determinant lemma for an invertible square matrix H where   \begin{eqnarray} \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) = \det ({\bf H})(1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}). \end{eqnarray} (17)The update uses the rank-one perturbation vectors u and v, selecting by row and column, to manipulate a single element of the matrix, a node of our mesh, by a small perturbation. A first-order approximation of the derivative of a determinant via a rank-one perturbation is then   \begin{eqnarray} \frac{\partial \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) }{\partial x_n} = \det ({\bf H})\, {\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}. \end{eqnarray} (18)To evaluate the change in our shape function’s determinant, we use the partial derivative of an absolute function ∂|H| = H∂H/|H| and the inverse determinant equivalence $$\det ({\bf H}^{-1})=\det ({\bf H})^{-1}$$ so that   \begin{eqnarray} \frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}= \frac{\partial \, |\det ({\mathbf {E}})^{-1}|}{\partial x_n} = \frac{\partial \, |\det ({\mathbf {E}}^{-1})|}{\partial x_n} \end{eqnarray} (19)  \begin{eqnarray} \phantom{\frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}}= \frac{\det ({\mathbf {E}}^{-1}) \, \frac{\partial (\det {\mathbf {E}}^{-1})}{\partial x_n} }{|\det ({\mathbf {E}}^{-1})|} \end{eqnarray} (20)and the partial derivative of the determinant   \begin{eqnarray} \frac{\partial \det ({\mathbf {E}}^{-1}) }{\partial x_n} = \det ({\mathbf {E}}^{-1})\, {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}} \end{eqnarray} (21)can be applied to (20) after reducing the determinants   \begin{eqnarray} \frac{\det ({\mathbf {E}}^{-1}) \det ({\mathbf {E}}^{-1})}{|\det ({\mathbf {E}}^{-1})|} = \frac{|\det {\mathbf {E}}|}{\det ({\mathbf {E}})^2} = \frac{1}{|\det {\mathbf {E}}|} \end{eqnarray} (22)so that   \begin{eqnarray} \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} &= \frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \end{eqnarray} (23) The partial derivative of the reduced shape matrix can be approximated using the Sherman–Morrison formula   \begin{eqnarray} ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} )^{-1} = {\bf H}^{-1} - \frac{{\bf H}^{-1} {\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}}{1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1} {\mathbf {u}}} \end{eqnarray} (24)to get a rank-one update   \begin{eqnarray} \frac{\partial {\mathbf {E}}_{\setminus 1} }{\partial x_n} = -({\mathbf {E}}{\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\mathbf {E}})_{\setminus 1} \end{eqnarray} (25)for a small perturbation such that $${\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}\ll 1$$. To go from a 2-D solution to a 2.5-D solution, a correction term k2T appears in the system matrices so that the shape matrices S(e) are extended to become S(e) + k2T(e) for spatial wave-number k.   \begin{eqnarray} {\mathbf {T}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (26)For the 2.5-D position Jacobian, this change adds a new partial derivative term ∂T(e)/∂xn  \begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} \, \rightarrow \, \frac{\partial ({\mathbf {S}}^{(e)}+k^2{\mathbf {T}}^{(e)})}{\partial x_n} = \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} + k^2 \frac{\partial {\mathbf {T}}^{(e)}}{\partial x_n} \end{eqnarray} (27)where ∂S(e)/∂xn is already available from the 2-D calculations, and the additional term ∂T(e)/∂xn may be derived   \begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (28)but we already have $$\partial |\det {\mathbf {E}}|^{-1} / \partial x_n$$ from deriving the partial derivatives of the S(e) term (23) giving   \begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (29)for a linear interpolatory shape function E perturbing a node at row u and column v. This adjoint or rank-one perturbation method for the 2.5-D electrode position Jacobian may be calculated much more quickly than a direct perturbation method because the system matrices do not need to be recalculated and inverted to determine the change in measurements due to electrode movement. The 2.5-D electrode position Jacobian (13) was found to be 25.9 times faster than the equivalent 3-D rank-one update method (Gómez-Laberge & Adler 2008), implementing (12) in 3-D, for mesh geometries used in this work (Intel Core i5-2500K 4-core processor at 3.30 GHz with 32 GB memory). 2-D electrode position Jacobian estimates differ significantly from 3-D solutions (Fig. 8), so have not been presented in Fig. 9. The computational cost of the 2-D rank-one update method for calculating the electrode position Jacobian (Gómez-Laberge & Adler 2008) was orders of magnitude faster than a naïve 2-D perturbation method. The 2-D rank-one update was 7.1 times faster than the 2.5-D method across most mesh sizes, which is accounted for by the numerical integration implied by (13). There are likely to be further gains from optimizing this implementation for multiple processing cores because key portions of the Jacobian calculation (14), (15), (26), (29) can be performed in parallel and the Jacobian typically consumes a significant portion of the total calculation time in each Gauss–Newton iteration (Boyle et al. 2012b). Mesh density was measured as the inverse of the maximum element height h (elements per metre) for both 2-D and 3-D meshes. The relative speed-up for a particular mesh density 1/h grows as a function of the number of FEM mesh nodes n and elements which must be calculated in the Jacobian where n = O(h−2) in 2-D and n = O(h−3) in 3-D; the larger the mesh the greater the benefit conferred by the 2.5-D Jacobian approach. Figure 9. View largeDownload slide The 2.5-D Jacobian calculations scale with FEM node count: the 2.5-D movement Jacobian speed advantage over a rank-one 3-D calculation (Gómez-Laberge & Adler 2008) grows as mesh density (FEM elements per metre) 1/h, where h is the maximum element height for the entire mesh. Error bars show max/min run times over 20 runs, the run time is closely related to the number of nodes in the FEM mesh where 2-D meshes have far fewer nodes for the same mesh density. Figure 9. View largeDownload slide The 2.5-D Jacobian calculations scale with FEM node count: the 2.5-D movement Jacobian speed advantage over a rank-one 3-D calculation (Gómez-Laberge & Adler 2008) grows as mesh density (FEM elements per metre) 1/h, where h is the maximum element height for the entire mesh. Error bars show max/min run times over 20 runs, the run time is closely related to the number of nodes in the FEM mesh where 2-D meshes have far fewer nodes for the same mesh density. 6 RESULTS The column ℓ2-norm sensitivity (the diagonal of $${{\bf J}}^{{\sf T}} {{\bf J}}$$) was plotted by replacing reconstructed resistivity with the log of estimated sensitivity. Sensitivity in these plots was expected to be greatest near the electrodes and diminish elsewhere. Simple Dirichlet boundary conditions away from the electrodes, used in these simulations, introduced errors, which were observable as variations in sensitivity at unexpected locations. We wished to model an approximately half-space forward model but unexpected increases in sensitivity near the sides and bottom were found to be caused by the boundary conditions which were deflecting current flow. Boundary condition errors can be corrected in a number of ways, the simplest of which is to increase the modelled domain until the error is small enough. One could, alternatively, implement appropriate ‘infinite elements’ at the boundary (Babus̆ka 1972). Another method is to estimate a ‘primary’ field for each stimulus using an analytic half-space model or very detailed one-time-use FEM mesh, and then calculate a ‘secondary’ field update on a smaller FEM with different resistivity as a correction (Günther et al. 2006). The primary-secondary type methods rely on small changes in resistivity far enough from the boundary to leave the ‘primary’ field largely unperturbed and it is not immediately obvious how this method may be affected by electrode movements or surface deformation without recalculating the primary field at each update. Neumann and mixed boundary conditions away from the electrodes were not considered. We have used an expanded model, in the interests of reliable results under deformed boundaries, at the expense of some lost computational efficiencies. Regions of high sensitivity were initially noted at depth where little sensitivity was expected. To determine how far the FEM model boundaries needed to be extended, an analytic PEM half-space model was compared to CEM homogeneous resistivity FEM simulated measurements. The model boundaries were extended approximately one electrode array length in each of the +x, −x and −z directions. The boundary extension reduced boundary condition related errors in simulated measurements to within measured noise levels and removed the artefacts from the sensitivity plots. The resistivity sensitivity S was plotted relative to the maximum sensitivity,   \begin{eqnarray} {\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V} \end{eqnarray} (30)using the conductivity Jacobian Jσ, measurement inverse covariance/weighting W, and element volumes $$\bf V$$, for line#1 (March 2008) and line#5 (February 2013) using the surveyed locations (Fig. 10). The region near the electrodes has been presented with annotations matching Fig. 6, as well as an image of the complete model. The initial and final resistivities were independently reconstructed using the surveyed locations (Fig. 11) and the difference between initial and final was used to create the expected resistivity change (Figs 12e and f). The initial resistivity for line#1 closely matched those published in Wilkinson et al. (2010) (Fig. 6) and achieved a similar <1 per cent RMS measurement misfit relative to a homogeneous resistivity model. Qualitatively and quantitatively, the two reconstructions (Figs 6 and 11a) are very similar. Figure 10. View largeDownload slide Inverse model parametrization with colours showing relative sensitivity $${\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V}$$ as S/max (S) for homogeneous resistivity; (a) line#1 March 2008, and (b) line#5 February 2013; note the distinct difference in slope profile between (a) and (b). Figure 10. View largeDownload slide Inverse model parametrization with colours showing relative sensitivity $${\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V}$$ as S/max (S) for homogeneous resistivity; (a) line#1 March 2008, and (b) line#5 February 2013; note the distinct difference in slope profile between (a) and (b). Figure 11. View largeDownload slide Resistivity-only reconstructions using true electrode positions as measured by RTK GPS; (a) line#1 March 2008 to (b) April 2009 (λ = 52.1, σ0 = 31.3 Ωm), and (c) line#5 February 2013 to (d) February 2014 (λ = 42.7, σ0 = 25.6 Ωm). Figure 11. View largeDownload slide Resistivity-only reconstructions using true electrode positions as measured by RTK GPS; (a) line#1 March 2008 to (b) April 2009 (λ = 52.1, σ0 = 31.3 Ωm), and (c) line#5 February 2013 to (d) February 2014 (λ = 42.7, σ0 = 25.6 Ωm). Figure 12. View largeDownload slide Change in resistivity and electrode movement for joint movement reconstructions (λσ = 0.1, λx = 0.07), (a–d) line#1, March 2008 to April 2009, and (e–h) line#5, February 2013 to February 2014. Figure 12. View largeDownload slide Change in resistivity and electrode movement for joint movement reconstructions (λσ = 0.1, λx = 0.07), (a–d) line#1, March 2008 to April 2009, and (e–h) line#5, February 2013 to February 2014. Electrode movements were initially reconstructed without allowing resistivity change using an independent implementation of the electrode movement iterative solver. The observed behaviour of the algorithm (Fig. 13) was to first minimize the error in the electrode spacing (corresponding to the largest measurement misfit), and then to approach a more ‘correct’ solution. Artefacts of this approach to the minimum over the optimization surface remain where large steps in the true electrode movement are reconstructed as a balanced step by adjacent electrodes that come close to the true displacement between electrodes near that movement. The joint inversion code was checked against this result by setting the movement-resistivity balance parameter β to strongly favour electrode movement. Reconstructions showed no resistivity change and movements that were in close agreement with the movement-only reconstruction code. Small variations still existed between the two results due to differences in the Gauss–Newton implementation and inexact line search. These variations were small with respect to the overall electrode movement solution. When electrode movements were reconstructed with resistivity changes (Fig. 12), some portion of the reconstructed electrode movement was lost in favour of reconstructed resistivity change. Figure 13. View largeDownload slide Electrode movement without allowing for resistivity changes, iterations for (a) line#1, March 2008 to April 2009, and (b) line#5 movement, February 2013 to February 2014. Figure 13. View largeDownload slide Electrode movement without allowing for resistivity changes, iterations for (a) line#1, March 2008 to April 2009, and (b) line#5 movement, February 2013 to February 2014. Due to the large electrode movements, it was found helpful to perform a crude version of successive relaxation. The first three iterations of the Gauss–Newton reconstruction were performed with an electrode movement hyperparameter that was an order of magnitude larger than following iterations. Without this adaptation, the reconstructed electrode movements showed poor agreement with measured locations, presumably because the Gauss–Newton iterations were trapped in a local minimum which favoured constructing resistivity artefacts near the electrodes. Exploring the space of hyperparameters near the selected hyperparameter did not reveal one which achieved better electrode movement reconstruction. 7 DISCUSSION Resistivity was reconstructed for measured initial and final electrode locations (Fig. 11) which serve as an ‘ideal’ reconstruction. Resistivity and electrode displacement were simultaneously reconstructed for a survey located on a slowly moving landslide. Results exhibit some measure of oscillatory artefacts in the reconstructed movement (Figs 12c and g). The resistivity distribution for line#1 (March 2008–April 2009; Figs 11a and b) changed by a relatively small amount when using true electrode locations before and after movement. This would suggest that, beyond the ground motion at the surface, no structural changes in the near surface seem to have occurred. It seems plausible that the increased area of low resistivity WMF might be indicative of increased saturation of the soil which led to the translational slide of WMF material moving over SSF substrate at the surface. The resistivity changes for line#5 (February 2013–February 2014; Figs 11c and d), show a significantly different distribution after ground movement, which is interesting given that the two electrode lines are within 40 m of each other. The line#5 measurements (Fig. 11d) occurred after a very wet summer and winter period where there was a lot of seepage at the base of the lobe causing the deeper reductions in resistivity between z = 40 and 60 m. At the surface of the lobes, resistivity increased due to cracking of the top layer. A cracked surface experienced accelerated evaporation over increased surface area, resulting in localized resistivity increase. To a lesser extent, areas affected by surface cracking also showed increased resistivity due to the change in topography: relatively little current would be conducted across the air gap in cracks, contributing to an average increase in bulk conductivity, but the cracks do increase the surface along which current flows resulting in an effective increase in resistivity. We take particular note of the change in SSF resistivity around x = 80 m which may have developed vertical connectivity between the overlying WMF and RMF below, allowing vertical drainage. The flow would be from the saturated WMF, along the WMF-SSF boundary to the surface, then to a region of vertical connectivity downward through the SSF (x = 80 m), and then into the RMF where it has pooled underground. This proposed flow path might also explain the increase in resistivity in the SSF at x = 60 m: if the vertical connectivity were in a roughly vertical plane, it would cut off the outer section of the SSF and that outer section would drain downwards into the RMF leading to an increase in resistivity. The deeper segment of the SSF (x > 80 m) would maintain its resistivity because the general connectivity and saturation have not changed by much. The resistivity change may also be induced by model error in electrode placement or topography: the 2.5-D model limits fidelity in some respects. The poor quality of the line#5 post-movement data, as measured by reciprocal error, may be the cause of these changes, though the locations of the reciprocal errors along the electrode array were distributed along the length of the array so that we expect no concentrated region of low sensitivity that could cause resistivity changes in the reconstruction such as those observed Fig. 11(d). Given this speculation, it would be interesting to investigate this potential vertical fault in the SSF. Perhaps it is an indication of a major ground movement, still to come, as the lower slope drainage has changed significantly. The change in drainage may also help to stabilize the lower slope by providing a drainage path at-depth which will allow surface material in the SSF to consolidate. This stabilizing effect has been observed on other nearby slopes in previous years. It has been postulated in Uhlemann et al. (2017) that reactivation of the slope at line#1 was stabilized due to slope movement which caused preferential flow paths to open, lowering pore pressures on the slip surface of the lobe, thereby stabilizing the lobe. For a half-space model with a homogeneous resistivity, electrode positions are not unique. A translation of the entire set of electrodes will give identical measurements. Likewise, a scaling of all electrode positions is equivalent to a change in the homogeneous resistivity. When conductivities are inhomogeneous the electrode locations are somewhat fixed by the locations of the inhomogeneities. Examples of electrode position non-uniqueness manifested itself in this data as large oscillations in the reconstructed electrode movement when no measures were taken to address the issue. Fixing the location of three electrodes at the upslope and downslope ends of the electrode array (Figs 12c and g) nearly eliminated these oscillations. As seen in the line#5 data, this is not necessarily a correct assumption, as both the top and bottom of a landslide may move, leading to resistivity artefacts. We infer that fixing these electrode locations was sufficient to damp the reconstructed movement’s oscillations because it fixes the relationship between a stimulus current, a measured voltage, and electrode separation (distance). Smoothing-type regularization of resistivity, used in this work, then controls how strictly the selected scaling of electrode separation is enforced. For example, electrodes that have contracted together in a region might cause a conductive artefact to be reconstructed under those electrodes. Increasing the resistivity regularization may suppress these artefacts and cause movement to be reconstructed by contracting the local electrode spacing to account for smaller than expected voltage measurements in the region. Both unscaled and log scaled electrode movement were tested and found to give results with similar absolute positional error. We have elected to present the unscaled electrode movement in our reconstructions as it is a less restrictive choice. One could imagine low angle slopes, wetlands or floodplains where the expected direction of movement would not be known a priori. Peat wetlands, muskeg, and most ground exposed to deep frost experience seasonal expansion and contraction, due to freeze-thaw cycles, water table changes, or water and gas accumulation and evaporation, which can result in uplift and ground shifting in directions other than downslope (Taber 1930; Hansell et al. 1983; Price 2003; Strack et al. 2006; Uhlemann et al. 2016c). A further reason to avoid dependence on the log movement scaling is the extension of this work to transverse movements where the restriction to movements only to one side of the array seems inappropriate. In the data sets examined here, there are transverse movements which were caused by material accumulating at the toe of the landslide and towards the edges of the earth flow. These transverse movements can be predicted for this particular data set based on the pre-existing topology: line#1 electrodes moved east, downhill into a gully, while the line#5 electrodes moved west, downhill into the same gully. Reconstructions for line#1 generally matched the true electrode locations within 0.20 m for movements of up to 1.46 m with the exception of electrode #9 and the three electrodes #6–#8 at the step in electrode position (Fig. 12a). Compared to Wilkinson et al. (2010) (≤0.2 m position error), these results are marginally less accurate. It seems probable that our results for the line#1 data differ from those of Wilkinson et al. (2010) due to the simultaneous resistivity and electrode position reconstruction, presented here, which removes artificial ordering constraints that occur with the sequential method of Wilkinson et al. (2010). Another source of differences in our results with respect to Wilkinson et al. (2010) is the restriction to downslope movements by a log parametrization in Wilkinson et al. (2010). As mentioned previously, trials of this log parametrization method in our simultaneous resistivity and electrode position inversion did not lead to improved results. Our results for the line#5 data appear to closely correlate with Wilkinson et al. (2016), where reconstructed electrode locations were generally within 0.2 m excepting some electrodes with errors up to 1.0 m, though results are not presented in as much detail in that case. In contrast to Wilkinson et al. (2016), our reconstructions do not address movements transverse to the electrode line. Movement reconstructions for line #5 do not appear to be particularly accurate, perhaps due to the more significant resistivity changes inferred in the reconstruction and more widespread translational failure of the slope which shifted electrodes over most of the resistivity structure (Fig. 12g). These might be addressed by identifying the covariance between movement and resistivity change within the joint reconstruction algorithm regularization. It is also possible that some of the error in reconstructed electrode position may be due to the FEM discretization. An approach such as the Fréchet method for electrode movement may help to alleviate such issues, though in general it produces the same solutions as a 3-D rank-one update method (Dardé et al. 2012; Boyle et al. 2017). Adjusting the relationship between resistivity and movement regularization β caused greater electrode displacement error as resistivity regularization was reduced. These movement magnitudes represented movement of up to 32 per cent of the average 4.73 m electrode spacing, exceeding the joint resistivity-movement methods of Soleimani et al. (2006) which was limited to movements of approximately 1 per cent of electrode spacing. 8 CONCLUSIONS This work demonstrated the practical application of a joint electrode movement-resistivity reconstruction using an iterative Gauss–Newton regularized solver. The electrode position Jacobian was calculated on the current resistivity at each iteration. Reconstructions show reasonable agreement with RTK GPS measured electrode locations, available resistivity estimates and geological structure. The initial reconstructed resistivity model, used as a starting point for the electrode movement and resistivity reconstruction, was in close agreement with prior work (Figs 6 and 11a) (Wilkinson et al. 2010). Reconstructed changes in resistivity (Fig. 12) showed considerable variation, particularly around the region at the toe of the landslide. These changes in resistivity could be indicative of water saturation changes due to water seepage at the toe of the landslide or other geological causes. Another possibility is that the resistivity changes represent artefacts due to transverse and normal components of electrode movement which were not accounted for in these reconstructions. We note that, in general, even when electrode displacements were not entirely accurate compared to true electrode positions, the error in the estimated and true distances between electrode positions was quite accurate after one or two Gauss–Newton iterations. Errors in electrode spacing were distributed fairly evenly across the electrode array, after which the displacements shifted towards their true positions in most cases. This suggests that a parametrization for electrode movement that encompasses electrode spacing may lead to improved outcomes. The results suggest that electrode grids are effective not only for resistivity monitoring, but also as a means of ground motion detection which may provide a cost-effective approach for landslide monitoring. We are encouraged by these results and expect that with new protocols which measure between electrode lines, higher quality electrode position reconstructions will be possible which incorporate normal, lateral, and transverse electrode movements, as observed in the data sets presented here. Acknowledgements Identification of the Hollin Hill site, geology, field work and electrical measurements were carried out by the British Geological Survey. Thanks go out to the researchers at the British Geological Survey, Geophysical Tomography Team for generously sharing their time and data. They were supported by the Natural Environment Research Council (NERC) and their contributions to this work are published with the permission of the Executive Director of the British Geological Survey. In fond memory of Steve Gibson, and with thanks to Josie Gibson. Code for the algorithms presented in this work extend EIDORS, an open source Matlab toolkit for impedance imaging inverse problems (Adler & Lionheart 2006; Adler et al. 2015, 2017), and use NetGen for meshing (Schöberl 1997). This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC). Footnotes 1 Leica System 1200 RTK GPS in ‘kinematic mode’ real-time correction achieves as much as 10 mm (Root Mean Squared (RMS)) horizontal and 20 mm vertical accuracies (Merritt 2014). 2 Longitudinal movement being movement inline with the electrodes and along the surface. 3 Matlab interp1(X,Y,Xq,‘pchip’) 1-D interpolant based on Hermite derivatives. REFERENCES Adler A., Lionheart W.R.B., 2006. Uses and abuses of EIDORS: an extensible software base for EIT, Physiol. Meas. , 27( 5), S25– S42. https://doi.org/10.1088/0967-3334/27/5/S03 Google Scholar CrossRef Search ADS PubMed  Adler A., Guardo R., Berthiaume Y., 1996. Impedance imaging of lung ventilation: do we need to account for chest expansion?, IEEE Trans. Biomed. Eng. , 43( 4), 414– 420. https://doi.org/10.1109/10.486261 Google Scholar CrossRef Search ADS PubMed  Adler A., Gaburro R., Lionheart W. R.B., 2011. Electrical impedance tomography, in Mathematical Methods in Imaging , chap. 14, pp. 601– 654, ed. Scherzer O., Springer Science+Business Media. Adler A., Boyle A., Crabb M.G., Gagnon H., Grychtol B., Lesparre N., Lionheart W.R.B., 2015. EIDORS Version 3.8, in International Conference on Biomedical Applications of Electrical Impedance Tomography , Neuchâtel, Switzerland. Adler A., Boyle A., Braun F., Crabb M.G., Grychtol B., Lionheart W. R.B., Tregidgo H.F.J., Yerworth R., 2017. EIDORS 3.9, in 18th International Conference on Biomedical Applications of Electrical Impedance Tomography , Dartmouth, USA. Babus̆ka I., 1972. The finite element method for infinite domains. I, Math. Comput. , 26( 117), 1– 11. https://doi.org/10.1090/S0025-5718-1972-0298969-2 Google Scholar CrossRef Search ADS   Boyle A., 2010. The effect of boundary shape deformation on two-dimensional electrical impedance tomography, Master’s thesis , Carleton University, Ottawa, Canada. Boyle A., 2016. Geophysical applications of electrical impedance tomography, PhD thesis , Carleton University, Ottawa, Canada. Boyle A., Adler A., 2010. Electrode models under shape deformation in electrical impedance tomography, in 11th Conf. Electrical Impedance Tomography , University of Florida, Gainesville, USA. Google Scholar CrossRef Search ADS   Boyle A., Adler A., 2011. Impact of electrode area, contact impedance and boundary shape on EIT images, Physiol. Meas. , 32( 7), 745– 754. https://doi.org/10.1088/0967-3334/32/7/S02 Google Scholar CrossRef Search ADS PubMed  Boyle A., Adler A., 2016. Modelling with 2.5D approximations, in 17th Conference on Electrical Impedance Tomography , Stockholm, Sweden. Boyle A., Adler A., Lionheart W.R.B., 2012a. Shape deformation in two-dimensional electrical impedance tomography, IEEE Trans. Med. Imaging , 31( 12), 2185– 2193. https://doi.org/10.1109/TMI.2012.2204438 Google Scholar CrossRef Search ADS   Boyle A., Borsic A., Adler A., 2012b. Addressing the computational cost of large EIT solutions, Physiol. Meas. , 33( 5), 787– 800. https://doi.org/10.1088/0967-3334/33/5/787 Google Scholar CrossRef Search ADS   Boyle A., Wilkinson P., Chambers J., Lesparre N., Adler A., 2014. Slope stability monitoring through impedance imaging, in 15th Conf. Electrical Impedance Tomography , Gananoque, Canada. Boyle A., Crabb M.G., Jehl M., Lionheart W.R.B., Adler A., 2017. Methods for calculating the electrode movement Jacobian for impedance imaging, Physiol. Meas. , 38( 3), 555– 574. https://doi.org/10.1088/1361-6579/aa5b78 Google Scholar CrossRef Search ADS PubMed  Carlson R., Fritsch F., 1989. An algorithm for monotone piecewise bicubic interpolation, SIAM J. Numer. Anal. , 26( 1), 230– 238. https://doi.org/10.1137/0726013 Google Scholar CrossRef Search ADS   Chambers J. et al.  , 2011. Three-dimensional geophysical anatomy of an active landslide in Lias Group mudrocks, Cleveland Basin, UK, Geomorphology , 125( 4), 472– 484. https://doi.org/10.1016/j.geomorph.2010.09.017 Google Scholar CrossRef Search ADS   Cormen T., Leiserson C., Rivest R., Stein C., 1990. Introduction to Algorithms , 3rd edn, The MIT Press. Dardé J., Hakula H., Hyvönen N., Staboulis S., 2012. Fine-tuning electrode information in electrical impedance tomography, Inverse Probl. Imaging , 6( 3), 399– 421. https://doi.org/10.3934/ipi.2012.6.399 Google Scholar CrossRef Search ADS   Dardé J., Hyvönen N., Seppänen A., Staboulis S., 2013a. Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography, SIAM J. Imaging Sci. , 6( 1), 176– 198. https://doi.org/10.1137/120877301 Google Scholar CrossRef Search ADS   Dardé J., Hyvönen N., Seppänen A., Staboulis S., 2013b. Simultaneous recovery of admittivity and body shape in electrical impedance tomography: an experimental evaluation, Inverse Probl. , 29( 8), 1– 16. Google Scholar CrossRef Search ADS   Dey A., Morrison H., 1979. Resistivity modelling for arbitrarily shaped two-dimensional structures, Geophys. Prospect. , 27, 106– 136. https://doi.org/10.1111/j.1365-2478.1979.tb00961.x Google Scholar CrossRef Search ADS   Friedel S., Thielen A., Springman S., 2006. Investigation of a slope endangered by rainfall-induced landslides using 3D resistivity tomography and geotechnical testing, J. appl. Geophys. , 60( 2), 100– 114. https://doi.org/10.1016/j.jappgeo.2006.01.001 Google Scholar CrossRef Search ADS   Gómez-Laberge C., Adler A., 2008. Direct EIT Jacobian calculations for conductivity change and electrode movement, Physiol. Meas. , 29( 6), S89– S99. https://doi.org/10.1088/0967-3334/29/6/S08 Google Scholar CrossRef Search ADS PubMed  Günther T., Rücker C., Spitzer K., 2006. Three-dimensional modelling and inversion of DC resistivity data incorporating topography - II. Inversion, Geophys. J. Int. , 166( 2), 506– 517. https://doi.org/10.1111/j.1365-246X.2006.03011.x Google Scholar CrossRef Search ADS   Hansell R., Scott P., Staniforth R., Svoboda J., 1983. Permafrost development in the intertidal zone at Churchill, Manitoba: a possible mechanism for accelerated beach uplift, Arctic , 36( 2), 198– 203. https://doi.org/10.14430/arctic2263 Google Scholar CrossRef Search ADS   Hyvönen N., Seppänen A., Staboulis S., 2014. Optimizing electrode positions in electrical impedance tomography, SIAM J. Appl. Math. , 74( 6), 1831– 1851. https://doi.org/10.1137/140966174 Google Scholar CrossRef Search ADS   Jehl M., Avery J., Malone E., Holder D., Betcke T., 2015. Correcting electrode modelling errors in EIT on realistic 3D head models, Physiol. Meas. , 36( 12), 2423– 2442. https://doi.org/10.1088/0967-3334/36/12/2423 Google Scholar CrossRef Search ADS PubMed  Jomard H., Lebourg T., Binet S., Tric E., Hernandez M., 2007. Characterization of an internal slope movement structure by hydrogeophysical surveying, Terra Nova , 19( 1), 48– 57. https://doi.org/10.1111/j.1365-3121.2006.00712.x Google Scholar CrossRef Search ADS   Jongmans D., Garambois S., 2007. Geophysical investigation of landslides: a review, Bull. Soc. France , 178( 2), 101– 112. https://doi.org/10.2113/gssgfbull.178.2.101 Google Scholar CrossRef Search ADS   Kim J., Yi M., Supper R., Ottowitz D., 2014. Simultaneous inversion of resistivity structure and electrode locations in ERT, in 20th European Meeting of Environmental and Engineering Geophysics: Near Surface 2014 , pp. 560– 564, Athens, Greece. Kuras O., Pritchard J., Meldrum P., Chambers J., Wilkinson P., Ogilvy R., Wealthall G., 2009. Monitoring hyrdaulic processes with automated time-lapse electrical resistivity tomography (ALERT), C. R. Geosci. , 341( 10–11), 868– 885. https://doi.org/10.1016/j.crte.2009.07.010 Google Scholar CrossRef Search ADS   Lapenna V., Lorenzo P., Perrone A., Piscitelli S., Rizzo E., Sdao F., 2005. 2D electrical resistivity imaging of some complex landslides in Lucanian Apennine chain, southern Italy, Geophysics , 70( 3), B11– B18. https://doi.org/10.1190/1.1926571 Google Scholar CrossRef Search ADS   Lebourg T., Binet S., Tric E., Jomard H., Bedoui S.E., 2005. Geophysical survey to estimate the 3D sliding surface and the 4D evolution of the water pressure on part of a deep seated landslide, Terra Nova , 17( 5), 399– 406. https://doi.org/10.1111/j.1365-3121.2005.00623.x Google Scholar CrossRef Search ADS   Lebourg T., Hernandez M., Zerathe S., Bedoui S.E., Jomard H., Fresia B., 2010. Landslides triggered factors analysed by time lapse electrical survey and multidimensional statistical approach, Eng. Geol. , 114( 3–4), 238– 250. https://doi.org/10.1016/j.enggeo.2010.05.001 Google Scholar CrossRef Search ADS   Loke M., 2017. RES2DINV: Rapid 2-D Resistivity & IP Inversion using the Least-squares Method , Geotomosoft Solutions. Loke M., Wilkinson P., Chambers J., 2015. Rapid inversion of data from 2-D and from 3-D resistivity surveys with shifted electrodes, in Near Surface Geoscience 2015-21st European Meeting of Environmental and Engineering Geophysics , Turin, Italy. Loke M., Wilkinson P., Chambers J., 2016. 3-D resistivity inversion with electrodes displacements, in 25th International Geophysical Conference and Exhibition , pp. 806– 810, Adelaide, Australia. Google Scholar CrossRef Search ADS   Loke M., Wilkinson P., Chambers J., Meldrum P., 2017. Rapid inversion of data from 2D resistivity surveys with electrode displacements, Geophys. Prospect. , doi:10.1111/1365-2478.12522. Merritt A., 2014. 4D geophysical monitoring of hydrogeological precursors to landslide activation, PhD thesis , University of Leeds, Leeds, United Kingdom. Merritt A. et al.  , 2014. 3D ground model development for an active landslide in lias mudrocks using geophysical, remote sensing and geotechnical methods, Landslides , 11( 4), 537– 550. https://doi.org/10.1007/s10346-013-0409-1 Google Scholar CrossRef Search ADS   Naudetb V., Lazzari M., Perrone A., Loperte A., Piscitelli S., Lapenna V., 2008. Integrated geophysical and geomorphological approach to investigate the snowmelt-triggered landslide of Bosco Piccolo village (Basilicata, southern Italy), Eng. Geol. , 98( 3–4), 156– 167. https://doi.org/10.1016/j.enggeo.2008.02.008 Google Scholar CrossRef Search ADS   Nocedal J., Wright S., 1999. Numerical Optimization , Springer Science Business Media. Google Scholar CrossRef Search ADS   Oldenborger G., Routh P., Knoll M., 2005. Sensitivity of electrical resistivity tomography data to electrode position errors, Geophys. J. Int. , 163( 1), 1– 9. https://doi.org/10.1111/j.1365-246X.2005.02714.x Google Scholar CrossRef Search ADS   Perrone A., Iannuzzi A., Lapenna V., Lorenzo P., Piscitelli S., Rizzo E., Sdao F., 2004. High-resolution electrical imaging of the Varco d’Izzo earthflow (southern Italy), J. Appl. Geophys. , 56( 1), 17– 29. https://doi.org/10.1016/j.jappgeo.2004.03.004 Google Scholar CrossRef Search ADS   Piegaria E., Cataudella V., Maio R.D., Milano L., Nicodemi M., Soldovieri M., 2009. Electrical resistivity tomography and statistical analysis in landslide modelling: a conceptual approach, J. Appl. Geophys. , 68( 2), 151– 158. https://doi.org/10.1016/j.jappgeo.2008.10.014 Google Scholar CrossRef Search ADS   Price J., 2003. Role and character of seasonal peat soil deformation on the hydrology of undisturbed and cutover peatlands, Water Resour. Res. , 39( 9), SBH–1– SBH– 10. Rücker C., Günther T., 2011. The simulation of finite ERT electrodes using the complete electrode model, Geophysics , 76( 4), F227– F238. https://doi.org/10.1190/1.3581356 Google Scholar CrossRef Search ADS   Sass O., Bell R., Glade T., 2008. Comparison of GPR, 2D-resistivity and traditional techniques for the subsurface exploration of the Öschingen landslide, Swabian Alb (Germany), Geomorphology , 93( 1-2), 89– 103. https://doi.org/10.1016/j.geomorph.2006.12.019 Google Scholar CrossRef Search ADS   Schöberl J., 1997. NETGEN an advancing front 2D/3D-mesh generator based on abstract rules, Comput. Vis. Sci. , 1( 1), 41– 52. https://doi.org/10.1007/s007910050004 Google Scholar CrossRef Search ADS   Soleimani M., Gómez-Laberge C., Adler A., 2006. Imaging of conductivity changes and electrode movement in EIT, Physiol. Meas. , 27( 5), S103– S113. https://doi.org/10.1088/0967-3334/27/5/S09 Google Scholar CrossRef Search ADS PubMed  Somersalo E., Cheney M., Isaacson D., 1992. Existence and uniqueness for electrode models for electric current computed tomography, SIAM J. Appl. Math. , 52( 4), 1023– 1040. https://doi.org/10.1137/0152060 Google Scholar CrossRef Search ADS   Strack M., Kellner E., Waddington J., 2006. Effect of entrapped gas on peatland surface level fluctuations, Hydrol. Process. , 20( 17), 3611– 3622. https://doi.org/10.1002/hyp.6518 Google Scholar CrossRef Search ADS   Supper R. et al.  , 2014. Geoelectrical monitoring: an innovative method to supplement landslide surveillance and early warning, Near Surf. Geophys. , 12( 1), 133– 150. Suzuki K., Higashi S., 2001. Groundwater flow after heavy rain in landslide-slope area from 2-D inversion of resistivity monitoring data, Geophysics , 66( 3), 733– 743. https://doi.org/10.1190/1.1444963 Google Scholar CrossRef Search ADS   Taber S., 1930. The mechanics of frost heaving, J. Geol. , 38( 4), 303– 317. https://doi.org/10.1086/623720 Google Scholar CrossRef Search ADS   Uhlemann S., Hagedorn S., Dashwood B., Maurer H., Gunn D., Dijkstra T., Chambers J., 2016a. Landslide characterization using P- and S-wave seismic refraction tomography — the importance of elastic moduli, J. Appl. Geophys. , 134( 1), 64– 76. https://doi.org/10.1016/j.jappgeo.2016.08.014 Google Scholar CrossRef Search ADS   Uhlemann S. et al.  , 2016b. Assessment of ground-based monitoring techniques applied to landslide investigations, Geomorphology , 253( 1), 438– 451. https://doi.org/10.1016/j.geomorph.2015.10.027 Google Scholar CrossRef Search ADS   Uhlemann S., Sorensen J., House A., Wilkinson P., Roberts C., Gooddy D., Binley A., Chambers J., 2016c. Integrated time-lapse geoelectrical imaging of wetland hydrological processes, Water Resour. Res. , 52( 3), 1607– 1625. https://doi.org/10.1002/2015WR017932 Google Scholar CrossRef Search ADS   Uhlemann S. et al.  , 2017. 4D imaging of moisture dynamics during landslide reactivation, J. geophys. Res. , 122( 1), 398– 418. https://doi.org/10.1002/2016JF003983 Google Scholar CrossRef Search ADS   Wagner F.M., Bergmann P., Rücker C., Wiese B., Labitzke T., Schmidt-Hattenberger C., Maurer H., 2015. Impact and mitigation of borehole related effects in permanent crosshole resistivity imaging: an example from the Ketzin CO2 storage site, J. Appl. Geophys. , 123( 1), 102– 111. https://doi.org/10.1016/j.jappgeo.2015.10.005 Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Lelliot M., Wealthall G., Ogilvy R., 2008. Extreme sensitivity of crosshole electrical resistivity tomography measurements to geometric errors, Geophys. J. Int. , 173( 1), 49– 62. https://doi.org/10.1111/j.1365-246X.2008.03725.x Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Meldrum P., Gunn D., Ogilvy R., Kuras O., 2010. Predicting the movements of permanently installed electrodes on an active landslide using time-lapse geoelectrical resistivity data only, Geophys. J. Int. , 183( 2), 543– 556. https://doi.org/10.1111/j.1365-246X.2010.04760.x Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Uhlemann S., Meldrum P., Smith A., Dixon N., Loke M., 2016. Reconstruction of landslide movements by inversion of 4-D electrical resistivity tomography monitoring data, Geophys. Res. Lett. , 43( 3), 1166– 1174. https://doi.org/10.1002/2015GL067494 Google Scholar CrossRef Search ADS   Zhou B., Dhalin T., 2003. Properties and effects of measurement errors on 2D resistivity imaging surveying, Near Surf. Geophys. , 1( 3), 105– 117. https://doi.org/10.3997/1873-0604.2003001 © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Jointly reconstructing ground motion and resistivity for ERT-based slope stability monitoring

, Volume 212 (2) – Feb 1, 2018
16 pages

/lp/ou_press/jointly-reconstructing-ground-motion-and-resistivity-for-ert-based-vh9l0hyYbw
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx453
Publisher site
See Article on Publisher Site

### Abstract

Summary Electrical resistivity tomography (ERT) is increasingly being used to investigate unstable slopes and monitor the hydrogeological processes within. But movement of electrodes or incorrect placement of electrodes with respect to an assumed model can introduce significant resistivity artefacts into the reconstruction. In this work, we demonstrate a joint resistivity and electrode movement reconstruction algorithm within an iterative Gauss–Newton framework. We apply this to ERT monitoring data from an active slow-moving landslide in the UK. Results show fewer resistivity artefacts and suggest that electrode movement and resistivity can be reconstructed at the same time under certain conditions. A new 2.5-D formulation for the electrode position Jacobian is developed and is shown to give accurate numerical solutions when compared to the adjoint method on 3-D models. On large finite element meshes, the calculation time of the newly developed approach was also proven to be orders of magnitude faster than the 3-D adjoint method and addressed modelling errors in the 2-D perturbation and adjoint electrode position Jacobian. Hydrogeophysics, Electrical resistivity tomography (ERT), Inverse theory, Tomography 1 INTRODUCTION Electrical resistivity tomography (ERT) is a geophysical method where current is injected between pairs of electrodes attached to the surface of the earth while the voltage differences are measured between other electrode pairs. From this data, one may reconstruct a resistivity distribution for the near surface earth that best fits the available data. For ERT, relatively small electrode movements or boundary modelling errors can lead to reconstructions with resistivity artefacts so severe that the resulting image is difficult to interpret (Zhou & Dhalin 2003). Similar reconstruction artefacts have been observed for biomedical electrical impedance tomography (EIT) where the fundamental mathematical problem is identical to ERT (Adler et al. 1996, 2011; Boyle & Adler 2011). Intermittent electrode movement is expected during long-term monitoring of an active landslide site. One could re-survey the electrode locations after each movement, but this would be time consuming and expensive, particularly for remote locations. Changes in landslide topography are frequently accompanied by changes in the resistivity distribution of the slope: changes in soil properties, such as water saturation, affect slope stability and resistivity. A resistivity reconstruction that does not account for electrode movement, when movements have occurred, will tend to fit the data poorly. Attempts to force fit a resistivity distribution will, in most cases, lead to images with severe resistivity artefacts which, without further information, could be misinterpreted. Simultaneously adjusting electrode position and resistivity can allow for better data fit, reduced resistivity artefacts, and indications of electrode movement. The methods described here separate the confounding factors, electrode movement and resistivity changes, which have long hindered ERT landslide monitoring attempts. Previous work on electrode movement has largely focused on 2-D electrode position changes in the plane of the electrodes. In 2-D, changes in boundary shape that are conformal do not result in artefacts (Boyle et al. 2012a). Modelling in 3-D is important because it reflects the finite size of the electrodes and the resulting current flow. For 3-D reconstructions, conformal deformations are limited to rotation, scaling and translation which would normally be identified through an appropriate site survey. 3-D resistivity reconstructions can be prohibitive to calculate for a given level of fidelity and numerical convergence. The so-called 2.5-D method (Dey & Morrison 1979) combines a 2-D parametrization of the resistivity and electrode positions with a 3-D model of current flow in the ground. The method can enable significant reductions in the computational requirements for meshing and current density calculations while producing accurate results with respect to equivalent 3-D models. Electrode movement and resistivity have been reconstructed for biomedical problems where changes both electrode movement (±1 per cent of electrode spacing) and resistivity (±50 per cent) have been relatively mild, enabling Gauss–Newton difference solutions (Soleimani et al. 2006; Jehl et al. 2015), and for large deformations on cylindrical domains with surface-normal deformations (Dardé et al. 2013a,b; Hyvönen et al. 2014) in simulation and on tank data. For geophysics, the resistivity contrasts are generally strong, varying by orders of magnitude. These strong contrasts combined with large electrode movements, much greater than 1 per cent of electrode spacing, necessitate an iterative solution because the combined effects of the electrode movement give highly nonlinear interactions. Similar approaches are also currently being developed by other researchers (Kim et al. 2014; Wagner et al. 2015; Loke et al. 2015, 2016, 2017). Results on field data have been somewhat mixed and data dependent. One approach to account for electrode position modelling errors is to alternate between an electrode position update and resistivity update (Wilkinson et al. 2010, 2016). Such an orthogonal greedy optimization approach may be slow to converge, or even fail to converge, because descent over the optimization surface is restricted to movement parallel to the axes (Cormen et al. 1990). Ideally, both electrode movement and resistivity should be reconstructed simultaneously so that updates can descend in the globally optimal direction across the objective function. In this work, joint electrode position and resistivity reconstruction methods were applied, in 2.5-D, to two data sets (Figs 1 and 2, line#1 and line #5) from an active landslide in the UK. An efficient and accurate method for calculating the electrode position Jacobian on a 2.5-D resistivity model was developed. This work is motivated by developments in Boyle & Adler (2010); Boyle (2010) for 2-D electrode movement, and Boyle et al. (2014); Boyle (2016); Boyle & Adler (2016) where this data set and preliminary results were presented. A 2-D cross-sectional Finite Element Method (FEM) model of the local slope topology was built, with electrodes modelled using the Complete Electrode Model (CEM) (Somersalo et al. 1992; Rücker & Günther 2011). We reconstruct electrode movement and resistivity in an iterative Gauss–Newton algorithm, showing that under certain conditions both can be simultaneously reconstructed. Figure 1. View largeDownload slide Slope failures at Hollin Hill; (a) rotational failures near the top of the slope above line#5, June 2015, (b) earth flow at the toe of a landslide where line#5 runs through mid-slope with electrodes throughout the slipped material, June 2015, (c) satellite image of the hillside (2016), showing four landslide ‘lobes’, five lines of 32 electrodes as of February 2014, and base station location.[(c) Satellite imagery ©2016 DigitalGlobe, Getmapping plc, Infoterra Ltd Bluesky.] Figure 1. View largeDownload slide Slope failures at Hollin Hill; (a) rotational failures near the top of the slope above line#5, June 2015, (b) earth flow at the toe of a landslide where line#5 runs through mid-slope with electrodes throughout the slipped material, June 2015, (c) satellite image of the hillside (2016), showing four landslide ‘lobes’, five lines of 32 electrodes as of February 2014, and base station location.[(c) Satellite imagery ©2016 DigitalGlobe, Getmapping plc, Infoterra Ltd Bluesky.] Figure 2. View largeDownload slide Electrode locations; (a) electrode locations for five lines of thirty-two electrodes each, line#1 (blue) to line#5 (green) as of February 2014, (b) where electrodes were 10 mm × 170 mm spikes of stainless steel selected for its conductivity, low cost and corrosion resistance. Figure 2. View largeDownload slide Electrode locations; (a) electrode locations for five lines of thirty-two electrodes each, line#1 (blue) to line#5 (green) as of February 2014, (b) where electrodes were 10 mm × 170 mm spikes of stainless steel selected for its conductivity, low cost and corrosion resistance. 2 BACKGROUND Resistivity imaging has been used in geophysical investigations of the behaviour and precursors of landslides and failure surfaces (Perrone et al. 2004; Lapenna et al. 2005; Lebourg et al. 2005; Jongmans & Garambois 2007; Naudetb et al. 2008; Sass et al. 2008). The technique is attractive because resistivity is strongly dependent on water saturation, fracturing, clay content and weathering which are all key factors in slope stability (Piegaria et al. 2009). Slopes may be monitored over time to observe changes in these key parameters using automated systems to collect and analyse data on a daily basis (Kuras et al. 2009; Lebourg et al. 2010; Supper et al. 2014). Difference images may show immediate changes in water saturation (Suzuki & Higashi 2001; Friedel et al. 2006; Jomard et al. 2007) but are limited in their ability to perform long-term monitoring due to background resistivity changes and electrode movements. As illustrated in Fig. 3, small amounts of electrode movement may introduce significant artefacts (Zhou & Dhalin 2003; Oldenborger et al. 2005; Wilkinson et al. 2008). These artefacts may be reduced by accounting for changes in electrode position (Uhlemann et al. 2017). Figure 3. View largeDownload slide Electrode movement artefacts; simulated reconstructions with a conductive and insulating target (rectangular and square outlines, respectively), each with two electrodes (electrode #2 and #12 of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (a) 0 per cent, (b) 5 per cent and (c) 25 per cent of electrode spacing on a 2-D half-space reconstruction (40 dB Signal-to-Noise Ratio (SNR), λ = 0.01, Laplace regularization, Wenner measurement pattern). Single or well separated electrode location errors introduce characteristic ‘ringing’ artefacts that can overwhelm resistivity-based information. [Reproduced, in part, from Boyle et al. (2017), fig. 2., with author permission.] Figure 3. View largeDownload slide Electrode movement artefacts; simulated reconstructions with a conductive and insulating target (rectangular and square outlines, respectively), each with two electrodes (electrode #2 and #12 of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (a) 0 per cent, (b) 5 per cent and (c) 25 per cent of electrode spacing on a 2-D half-space reconstruction (40 dB Signal-to-Noise Ratio (SNR), λ = 0.01, Laplace regularization, Wenner measurement pattern). Single or well separated electrode location errors introduce characteristic ‘ringing’ artefacts that can overwhelm resistivity-based information. [Reproduced, in part, from Boyle et al. (2017), fig. 2., with author permission.] An active landslide was identified in North Yorkshire, UK (54°06΄39.2″N, 0°57΄34.9″W) and has been monitored since 2008 (Chambers et al. 2011; Wilkinson et al. 2016). The landslide is a very slow to slow moving composite multiple earth slide-earth flow (Uhlemann et al. 2016b), where a central portion of the slope has moved downhill by up to 3.5 m a year in some instances (Fig. 1). The slope itself exposes four formations: the Dogger Formation (DF), Whitby Mudstone Formation (WMF), Staithes Sandstone Formation (SSF) and Redcar Mudstone Formation (RMF), from top to bottom. The interfaces between sedimentary layers lie horizontally, with a gentle 5° dip to the North, determined through comparison of material interfaces at surrounding exposed slopes in the region (Chambers et al. 2011). The WMF, as the name implies, is a mudstone clay-based rock that is highly weathered and prone to movement during peak water saturation periods at Hollin Hill, occurring annually in the winter through early spring wet-season. The underlying SSF and unweathered RMF are structurally more competent, and landsliding is postulated to occur within the WMF (Uhlemann et al. 2016b). The slope lies at an average angle of 14° over a change in elevation of 40 m. The landslide movement is determined by translational movements upslope of the WMF–SSF interface. These lead to a loss of support of the local slope of the back-scarp, causing rotational failure. The mass accumulation at the WMF–SSF interface, driven by these translational movements, and elevated pore pressures then cause flow movements which form lobes (Uhlemann et al. 2016a,b). A grid of five rows of thirty-two permanently installed electrodes travelled along with these movements, shifting their positions relative to each other (Fig. 2). An automated impedance imaging survey was executed bi-daily and data were transmitted to a remote office for storage and analysis. In 2008–2009, a middle section of line#1 exhibited a translational failure with movements of up to 1.6 m. In 2013–2014, upper and middle segments of line#5 had rotational and translational failures of up to 3.5 m. The dipole-dipole measurement protocol used for line#1 and line#5 are visualized in Figs 4 and 5, showing the sequence of quadrupolar measurements with stimulus dipoles in red and measurement dipoles in blue. A single difference measurement is captured as one row of the protocol in the figure. In the adjacent vertical strip chart, the horizontally aligned measurements at the initial time (green) is contrasted with the homogeneous resistivity estimate (red) of what those measurements would be, and the change in measurements from initial to final measurement (blue). The rightmost strip chart shows the estimated error for each measurement as estimated from differences between reciprocal pairs of data for the initial measurements (green) and final measurements (blue). The figure serves to illustrate three challenges with this data set. First, the measurements do not fit, or nearly fit, a homogeneous resistivity model. Second, the change between the initial and final measurements (resistivity and movement) is as much as that caused by the inhomogeneous resistivity distribution in the initial measurements (resistivity only) for some data. Third, the measurement error varies over time so that the initial and final measurements have different associated error estimates. The sequence of measurements for line#5 (Fig. 5) started at the top of the slope and ran to the base, but is equivalent to that shown for line#1 (Fig. 4). The final measurements for line#5 have significantly worse reciprocal errors, but are the best of many data sets examined for a post-movement data set on line#5. Figure 4. View largeDownload slide Dipole–dipole measurement protocol for line#1; March 2008 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 32.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 4. View largeDownload slide Dipole–dipole measurement protocol for line#1; March 2008 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 32.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 5. View largeDownload slide Dipole–dipole measurement protocol for line#5; February 2013 and February 2014 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 26.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. Figure 5. View largeDownload slide Dipole–dipole measurement protocol for line#5; February 2013 and February 2014 measurements, (left) stimulus in red and measurements in blue, one row per difference measurement, (middle) initial difference measurements Va (green) compared to homogeneous resistivity at 26.1 Ωm shown as Va − Vh (red), and the change from initial to final measurements Vb − Va (blue), and (right) the reciprocal standard error as a percentage of the measurements estimated by comparing reciprocal measurements for initial (green) and final (blue) measurements. 3 PRIOR RECONSTRUCTIONS AT HOLLIN HILL Electrode movement has been previously reconstructed for line#1 data (2008/2009) along with separate resistivity sections (Wilkinson et al. 2010). The algorithm used in that instance achieved estimates of within 0.2 m (4 per cent of electrode spacing) of the electrode’s true positions as measured by Real-Time Kinematic Global Positioning System (RTK GPS).1 An initial resistivity-only reconstruction with known electrode positions gave a plausible distribution that was in good agreement with available geological evidence: borehole and auger data, evaluation of local geology, aerial LiDAR, differential GPS, lab correlation of representative samples to measured conductivities, piezometric pore pressure measurements, and in-situ rainfall and temperature records (Fig. 6; Chambers et al. 2011; Merritt et al. 2014). Resistivity reconstructions of data collected after movement exhibited artefacts. These artefacts were reduced when reconstructed movements were incorporated. The electrode movement was penalized in the upslope direction. The position Jacobian was estimated based on an analytic half-space model with a homogeneous resistivity assigned to each group of measurements based on electrode separation. The electrode movement was then reconstructed by minimizing   \begin{eqnarray} \arg \min \sqrt{\sum _i |e_i|^2} + \alpha \sum _j |m_j| + \beta \sum _j \theta (m_j)|m_j| \end{eqnarray} (1)for weighting factors α = 0.06 m−1, β = 0.32 m−1, Heaviside step function θ, movement mj, and misfit ei = rb − ra as the difference between predicted rb and observed ra apparent resistivity ratios. The apparent resistivity ratio was calculated as the ratio of the analytic half-space models (9) before and after movement   \begin{eqnarray} r = \frac{\rho ^{\prime }}{\rho } \left( \frac{ \frac{1}{AM^{\prime }} - \frac{1}{BM^{\prime }} - \frac{1}{AN^{\prime }} - \frac{1}{BN^{\prime }}}{\frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} - \frac{1}{BN}} \right) \end{eqnarray} (2)for homogeneous resistivity ρ, electrodes spaced AM, BM, AN, BN, and the updated locations and resistivity after movement indicated by primed ΄ variables. Dipole-dipole data for measurements n = 1 were discarded, as they were judged to be too dependent on transverse movements which were not reconstructed. A similar approach was applied (Wilkinson et al. 2016), to reconstruct 2-D surface xy-movements for the whole electrode grid by allowing for transverse movements through an additional weighting term   \begin{eqnarray} \arg \min &&\sum _i e_i^2 + \alpha \sum _j ||m_j|| + \beta \sum _j \theta (m_j^{(y)})|m_j^{(y)}| \nonumber\\ &&\quad+ \, \gamma \sum _j \theta (m_j^{(x)})|m_j^{(x)}| \end{eqnarray} (3)facilitating balancing the sensitivities of transverse $$m_j^{(y)}$$ and longitudinal $$m_j^{(x)}$$ movements by adjusting the weighting ratio β/γ. Figure 6. View largeDownload slide [From Wilkinson et al. (2010), figure 2, x, z-coordinates corrected:] 2-D resistivity image inverted from the baseline data set (2008 March). The inferred boundaries between the Whitby (WMF), Staithes (SSF) and Redcar (RMF) formations are shown by dotted black lines. Stratigraphic logs of boreholes are shown in grey scale. The main scarp and slipped WMF material are indicated by the black arrows. [Reproduced from Wilkinson et al. (2010), fig. 2, for comparison with Fig. 11(a) in this work. The RES2DINV reconstruction region, elsewhere in this paper referred to as the ‘RES2DINV outline,’ is selected by the software based on a heuristic pseudo-section method (Loke 2017).] Figure 6. View largeDownload slide [From Wilkinson et al. (2010), figure 2, x, z-coordinates corrected:] 2-D resistivity image inverted from the baseline data set (2008 March). The inferred boundaries between the Whitby (WMF), Staithes (SSF) and Redcar (RMF) formations are shown by dotted black lines. Stratigraphic logs of boreholes are shown in grey scale. The main scarp and slipped WMF material are indicated by the black arrows. [Reproduced from Wilkinson et al. (2010), fig. 2, for comparison with Fig. 11(a) in this work. The RES2DINV reconstruction region, elsewhere in this paper referred to as the ‘RES2DINV outline,’ is selected by the software based on a heuristic pseudo-section method (Loke 2017).] 4 METHOD While the work of Wilkinson et al. (2010, 2016), described in the previous section, focused on a sequential inversion procedure, here we simultaneously reconstruct resistivity and electrode position in a Gauss–Newton iterative framework. A 2-D model of the slope profile was constructed with independent parameters for resistivity and electrode position (Fig. 7). The ℓ2-norm of the data discrepancy and regularization terms was minimized by balancing the sensitivity of the resistivity parameters ρ and electrode movement x against the regularization terms   \begin{eqnarray} (\hat{\mathbf {\rho }},\hat{{\mathbf {x}}}) &=& \arg \min _{\mathbf {\rho },{\mathbf {x}}} \ ||\mathcal {F}(\mathbf {\rho },{\mathbf {x}}) - {\mathbf {m}}||_{\mathbf {W}} \nonumber \\ &&\quad + \, ||\lambda \beta {\mathbf {R}}_\mathbf {\rho }\log _{10}(\mathbf {\rho }_\star /\mathbf {\rho })||_{2} \nonumber \\ && \quad + \, || \lambda \eta {\mathbf {R}}_{{\mathbf {x}}} ({\mathbf {x}}- {\mathbf {x}}_\star )||_{2}. \end{eqnarray} (4)The optimal solution $$(\hat{\mathbf {\rho }},\hat{{\mathbf {x}}})$$ minimizes the data discrepancy between a forward model $$\mathcal {F}$$ and measurements m combined with some regularization Rρ and Rx acting to smooth an otherwise ill-posed and ill-conditioned inverse problem. The forward model was the parametrized 2.5-D model (Fig. 7) of resistivity ρ and longitudinal2 electrode position x and producing an estimate of expected measurements given the measurement protocol. Measurement error was weighted W based on estimates of measurement reliability. The regularization terms penalized changes in resistivity and electrode position from prior estimates (ρ⋆, x⋆). The relative sensitivity of the two types of parameters, resistivity and electrode movement, were balanced by adjusting the ratio of the scalar β/η. The overall strength of the regularization was adjusted by scaling both terms by λ. The resistivity was solved under a positive conductivity σ constraint by converting to inverse log units log 10σ = log 101/ρ. Figure 7. View largeDownload slide The 2.5-D forward model FEM parametrization; right-to-left expanding from the region surrounding a single electrode, to the scale of the electrode array, to the scale of the region surrounding the electrode array, the forward model is parametrized for electrode position x and resistivity ρ, resistive regions selected to demonstrate mesh structure. Figure 7. View largeDownload slide The 2.5-D forward model FEM parametrization; right-to-left expanding from the region surrounding a single electrode, to the scale of the electrode array, to the scale of the region surrounding the electrode array, the forward model is parametrized for electrode position x and resistivity ρ, resistive regions selected to demonstrate mesh structure. The objective function (4) was solved using the well-known iterative Gauss–Newton approach (Nocedal & Wright 1999,  section 10.2). The Gauss–Newton approach starts from an initial estimate (ρ0, x0), estimates the local slope of the objective function as the Jacobian (Jσ, Jx) to determine a new search direction (δρ, δx), and then performs an approximate line search in that direction to estimate an optimal step length α. At each iteration, the parameters were updated   \begin{eqnarray} (\delta \mathbf {\rho }, \delta {\mathbf {x}}) = -({{\bf J}}^{{\sf T}} {\mathbf {W}}{{\bf J}}+ \lambda ^2 {\bf Q})^{-1}({{\bf J}}^{{\sf T}} {\mathbf {W}}{\mathbf {b}}\, + \lambda ^2 {\bf Q} {\bf c}) \end{eqnarray} (5)  \begin{eqnarray} &&{(\mathbf {\rho }_{n+1}, {\mathbf {x}}_{n+1}) = (\mathbf {\rho }_n, {\mathbf {x}}_n) + \alpha (\delta \mathbf {\rho }, \delta {\mathbf {x}})} \\ &&{{\rm for }\, \,{\mathbf {b}}= \mathcal {F}(\mathbf {\rho }_n,{\mathbf {x}}_n) - {\mathbf {m}}}\nonumber \\ &&{{\bf c} = \left[ \begin{array}{c}\log _{10} 1/\mathbf {\rho }_n - \log _{10} 1/\mathbf {\rho }_\star \\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array} \right] = \left[ \begin{array}{c}\log _{10} \mathbf {\rho }_\star /\mathbf {\rho }_n\nonumber\\ {\mathbf {x}}_n-{\mathbf {x}}_\star \end{array}\right]}\nonumber \\ &&{{{\bf J}}= \left[ \begin{array}{c}\frac{-1}{\mathbf {\rho }_n \ln (10)} \, {{\bf J}}_{\sigma _n} \\ {{\bf J}}_{\mathbf {x}}\end{array} \right]} \nonumber \\ &&{{\mathbf {R}}= \left[ \begin{array}{c{@}{\quad}c} \beta {\mathbf {R}}_\mathbf {\rho } & 0 \\ 0 & \eta {\mathbf {R}}_{\mathbf {x}}\end{array}\right] \qquad {\rm and \, \, \, {\mathbf {R}}^{{\sf T}} {\mathbf {R}}= {\bf Q}}}\nonumber \end{eqnarray} (6)based on the data discrepancy b and the distance from the prior estimate $$\bf c$$ in combination with the regularization R, the slope of the objective function J and measurement weighting W. This formulation agrees with that of Boyle et al. (2017), but is modified to address a log resistivity parametrization. In contrast to a typical resistivity-only inversion, the reconstruction parameters, Jacobian and regularization R have been extended to incorporate the new electrode position parameters. The resistivity Jacobian has been calculated for conductivity $${{\bf J}}_{\sigma _n}$$ and then scaled, which is exactly equivalent to calculating the Jacobian on the log of resistivity. The movement Jacobian was found to be sensitive to resistivity changes in the small elements adjacent to the electrodes. To address this sensitivity, the log conductivity regularization combined a smoothing prior near electrodes with Tikhonov regularization away from the electrodes   \begin{eqnarray} {\mathbf {R}}_\rho = {\mathbf {I}}+ \nu {\mathbf {R}}_{\mathfrak {L}} \end{eqnarray} (7)  \begin{eqnarray} \nu = \exp -\, \frac{|x_e - x_\ell |}{|\bar{x}_{\ell \ell }|} \end{eqnarray} (8)where the Laplace smoothing $${\mathbf {R}}_{\mathfrak {L}}$$ was scaled ν by the distance between each FEM element centre xe and the nearest electrode xℓ, scaled by the average distance between electrodes $$\bar{x}_{\ell \ell }$$. This prior encourages small changes from the expected resistivity in regions with little information. In regions near the electrodes, changes will be pushed towards the spaces between electrodes rather than directly under the electrodes, as well as encouraging smooth transitions in resistivity near the electrodes. The regularization for movement was the Tikhonov prior Rx = I. In principle, there are correlated changes between resistivity and electrode movement (Kim et al. 2014) which may be partially accounted for by setting the off-diagonal blocks of the regularization matrix to non-zero values, but in practice these were not characterized and in the absence of a better guess were set to zero. The forward model was constructed as a 2-D cross-section based on the original electrode locations and the mesh was then perturbed by PCHIP3 interpolation (Carlson & Fritsch 1989) for electrode displacements. Forward modelled measurements and the Jacobians were calculated using the 2.5-D method (Dey & Morrison 1979). We have made use of the log conductivity to restrict the resistivity reconstruction to physically meaningful positive values. Experiments with the log movement constraint, to restrict electrodes to downslope movement, resulted in similar reconstructions to the ones presented here which used unrestricted electrode movement. For the log movement parametrization, the behaviour at each iteration was different to an unscaled movement due to the structure of the movements in this data set. Because each electrode that moved had a different magnitude of movement, the log movement reconstruction tended to solve for each electrode’s reconstructed displacement separately: one electrode per iteration. The apparent single electrode updates were actually an artefact of the log scaling, where smaller movements were reduced by orders of magnitude, so as to be inconsequential. Once the largest electrode placement error had been corrected, the next largest error would be addressed. We feel that this highlights the importance of careful selection of the reconstruction parametrization. It is possible that some Fourier decomposition or other basis of electrode movement with appropriate regularization might achieve greater reconstruction accuracy without artificially fixing any single electrode’s location or degrees of freedom for movement. 5 2.5-D POSITION JACOBIAN The 2-D electrode position Jacobian suffers from significant errors (Fig. 8), when compared to data from a 3-D model, which makes it inappropriate for 3-D problems. A 3-D electrode position Jacobian becomes prohibitively expensive to calculate as the mesh density grows. A 2.5-D approach offers a compromise by restricting sensitivity parametrization and electrode positions to the plane, while achieving high fidelity to equivalent 3-D models, at a fixed multiple of the 2-D computational effort. Figure 8. View largeDownload slide Jacobian sensitivity $${\rm diag}\,\,({{\bf J}}_x^{{\sf T}} {{\bf J}}_x)^\frac{1}{2}$$ on a 16-electrode, homogeneous (σ = 1) half-space model; 2-D rank-one electrode position Jacobian (Gómez-Laberge & Adler 2008) shows orders-of-magnitude error in sensitivity estimate, while 3-D analytic (9) and rank-one estimates (Gómez-Laberge & Adler 2008) are in close agreement with the 2.5-D (13) estimate. Figure 8. View largeDownload slide Jacobian sensitivity $${\rm diag}\,\,({{\bf J}}_x^{{\sf T}} {{\bf J}}_x)^\frac{1}{2}$$ on a 16-electrode, homogeneous (σ = 1) half-space model; 2-D rank-one electrode position Jacobian (Gómez-Laberge & Adler 2008) shows orders-of-magnitude error in sensitivity estimate, while 3-D analytic (9) and rank-one estimates (Gómez-Laberge & Adler 2008) are in close agreement with the 2.5-D (13) estimate. Two alternate methods were evaluated before developing the 2.5-D position Jacobian: a 2.5-D perturbation method which was relatively slow, and an analytic model of movement which was restricted to a homogeneous resistivity and the Point Electrode Model (PEM). Motivated by the short comings of these two methods, we develop the 2.5-D position Jacobian which is efficient and accounts for resistivity variation in the model using the CEM. The perturbation method used the underlying 2.5-D forward simulations of the full FEM resistivity model using the CEM. These movement perturbation calculations proved to be prohibitively slow because a mesh perturbation resulted in recalculations of the system matrices and a new inversion of that system matrix. The cost grows with number of electrodes and movement dimensions, so that for n = 32 electrodes, estimated in d = 1 dimensions, nd + 1 forward simulations were required at each iteration of the Gauss–Newton inverse solution. A line search typically required three to six forward simulations, meaning that the perturbation Jacobian required far more time to calculate than the rest of each iteration. When FEM mesh density was increased sufficiently to achieve good estimates of electrode position changes, the calculations took an unreasonable amount of time. Low mesh densities sped up the calculations but exhibited significant errors when compared to a 3-D perturbation solution at sufficient mesh density. An alternate solution was implemented by adapting the half-space analytic PEM forward model. For a half-space, the potential difference V measured over a homogeneous resistivity ρ with current I driven on the stimulus electrodes, is given by   \begin{eqnarray} V = \frac{\rho I}{2 \pi } \left( \frac{1}{AM} - \frac{1}{BM} - \frac{1}{AN} + \frac{1}{BN} \right) \end{eqnarray} (9)where each distance AM, BM, AN, BN is between a stimulus electrode and a measurement electrode. The model may be applied for any arbitrary pair-wise electrode placement. A position Jacobian may be constructed by applying the electrode movement perturbation. A next logical step would be to take the derivative of (9) with respect to electrode position and build up a solution specific approximate block-wise model of resistivity away from electrodes (Wilkinson et al. 2010), though this was not implemented in this work. Electrode positions were captured from the FEM model. A homogeneous resistivity was assigned based on the average resistivity of the current FEM model. The electrode position Jacobian produced by the half-space analytic perturbation method was compared to the 2.5-D perturbation Jacobian under homogeneous conditions. Using the modified half-space analytic perturbations was much faster to calculate than the 2.5-D perturbation method and reasonably accurate: the effect of topography was found to be somewhat accounted for. The loss of accuracy due to changing electrode models (CEM to PEM) and using a homogeneous resistivity were not so disruptive as to change signs in the position Jacobian although magnitudes were inaccurate. Motivated by the efficiency of the 2-D position Jacobian of Gómez-Laberge & Adler (2008), and the errors introduced by using a 2-D electrode position Jacobian, the 2.5-D position Jacobian was developed as the corollary of the 2.5-D conductivity Jacobian. In general, the 2.5-D forward solver is a well-known technique and is commonly used in geophysics ERT applications. An approximately half-space geometry, and a resistivity that is nearly uniform along one axis, fit well with a 2.5-D model, and occur naturally in many geological settings. By adapting the adjoint, or ‘standard method,’ of calculating the conductivity Jacobian by a rank-one update, to the 2.5-D technique, resistivity may be efficiently reconstructed. To reconstruct electrode movement, we desire a similar 2.5-D implementation for the electrode position Jacobian. In the following, we outline the 2.5-D conductivity Jacobian and present a new derivation for the 2.5-D electrode position Jacobian. We use the formulation of the 2-D conductivity Jacobian from Boyle et al. (2017) as a basis for these developments. The 2-D and 2.5-D conductivity Jacobians Jσ were calculated   \begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} {\mathbf {S}}\frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}{\mathbf {X}} \end{eqnarray} (10)  \begin{eqnarray} {{\bf J}}_{\sigma ,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} ({\mathbf {S}}+ k^2 {\mathbf {T}}) \frac{\partial {\mathbf {D}}}{\partial \sigma _n} {\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray} (11)for measurement selection $${\mathcal {T}}$$, system matrix A, mesh connectivity matrix C, mesh shape functions S, a conductivity change ∂D/∂σn, and the nodal voltages X = A−1B for stimulus B over an electrode modelled as a shunt in the y-direction when the 2-D FEM is meshed over the x–z plane. The system matrix $${\mathbf {A}}= {\mathbf {C}}^{{\sf T}} {\mathbf {S}}{\mathbf {D}}{\mathbf {C}}$$ is assembled from a connectivity matrix C mapping global node numbers to element-local node numbers, the element shape functions S and the conductivity D per element. For the 2.5-D position Jacobian, the system matrix Ak is specific to the spatial frequency k, as are the nodal voltages $${\mathbf {X}}_k = {\mathbf {A}}_k^{-1} {\mathbf {B}}$$. A perturbation node is selected at row u, column v, affecting a linear interpolatory shape function E. The 2.5-D position Jacobian may be calculated as an extension of the 2-D Jacobian, in an analogous way to (10), (11), as   \begin{eqnarray} {{\bf J}}_{x,{\rm 2D}} = - {\mathcal {T}}{\mathbf {A}}^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial {\mathbf {S}}}{\partial x_n} {\mathbf {D}}{\mathbf {C}}{\mathbf {X}} \end{eqnarray} (12)  \begin{eqnarray} {{\bf J}}_{x,{\rm 2.5D}} = - \frac{2}{\pi } \int _k {\mathcal {T}}{\mathbf {A}}_k^{-1} {\mathbf {C}}^{{\sf T}} \frac{\partial ({\mathbf {S}}+ k^2 {\mathbf {T}})}{\partial x_n} {\mathbf {D}}{\mathbf {C}}\frac{{\mathbf {X}}_k}{2} \end{eqnarray} (13)where the 2-D position Jacobian may be efficiently calculated using the rank-one update for the conductivity Jacobian (Gómez-Laberge & Adler 2008) with some new terms. Again, based on the 2-D formulation from Boyle et al. (2017), the element shape functions for element (e) may be summarized as the element shape matrix   \begin{eqnarray} {\mathbf {S}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|}{\mathbf {E}}_{\setminus 1}^{{\sf T}} {\mathbf {E}}_{\setminus 1} \end{eqnarray} (14)for a shape matrix E and a row-reduced version E∖1 where the top row of the matrix is removed. The shape matrix is distorted by having its nodes perturbed leading to the first-order estimate of element deformation   \begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} &=& \frac{1}{2}\, \Bigg ( \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} {\mathbf {E}}_{\setminus 1}^{{{\sf T}} } {\mathbf {E}}_{\setminus 1} \, \nonumber\\ &&\qquad + \frac{1}{|\det {\mathbf {E}}|} \bigg (\frac{\partial {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }}{\partial x_n} {\mathbf {E}}_{\setminus 1} + {\mathbf {E}}_{\setminus 1}^{{{\sf T}} }\frac{\partial {\mathbf {E}}_{\setminus 1}}{\partial x_n} \bigg ) \Bigg ) \end{eqnarray} (15)where xn refers to a global node numbered n and affects all elements e connected to that node. The local shape functions of each element, for first-order interpolatory shape functions on a 2-D mesh, are   (16)for a triangular element (blue) with three nodes p1, p2, p3. To calculate the partial derivatives of the first-order interpolatory shape functions, we make use of the matrix determinant lemma for an invertible square matrix H where   \begin{eqnarray} \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) = \det ({\bf H})(1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}). \end{eqnarray} (17)The update uses the rank-one perturbation vectors u and v, selecting by row and column, to manipulate a single element of the matrix, a node of our mesh, by a small perturbation. A first-order approximation of the derivative of a determinant via a rank-one perturbation is then   \begin{eqnarray} \frac{\partial \det ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} ) }{\partial x_n} = \det ({\bf H})\, {\mathbf {v}}^{{\sf T}} {\bf H}^{-1}{\mathbf {u}}. \end{eqnarray} (18)To evaluate the change in our shape function’s determinant, we use the partial derivative of an absolute function ∂|H| = H∂H/|H| and the inverse determinant equivalence $$\det ({\bf H}^{-1})=\det ({\bf H})^{-1}$$ so that   \begin{eqnarray} \frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}= \frac{\partial \, |\det ({\mathbf {E}})^{-1}|}{\partial x_n} = \frac{\partial \, |\det ({\mathbf {E}}^{-1})|}{\partial x_n} \end{eqnarray} (19)  \begin{eqnarray} \phantom{\frac{\partial \, |\det {\mathbf {E}}|^{-1} }{\partial x_n}}= \frac{\det ({\mathbf {E}}^{-1}) \, \frac{\partial (\det {\mathbf {E}}^{-1})}{\partial x_n} }{|\det ({\mathbf {E}}^{-1})|} \end{eqnarray} (20)and the partial derivative of the determinant   \begin{eqnarray} \frac{\partial \det ({\mathbf {E}}^{-1}) }{\partial x_n} = \det ({\mathbf {E}}^{-1})\, {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}} \end{eqnarray} (21)can be applied to (20) after reducing the determinants   \begin{eqnarray} \frac{\det ({\mathbf {E}}^{-1}) \det ({\mathbf {E}}^{-1})}{|\det ({\mathbf {E}}^{-1})|} = \frac{|\det {\mathbf {E}}|}{\det ({\mathbf {E}})^2} = \frac{1}{|\det {\mathbf {E}}|} \end{eqnarray} (22)so that   \begin{eqnarray} \frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} &= \frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \end{eqnarray} (23) The partial derivative of the reduced shape matrix can be approximated using the Sherman–Morrison formula   \begin{eqnarray} ({\bf H} + {\mathbf {u}}{\mathbf {v}}^{{\sf T}} )^{-1} = {\bf H}^{-1} - \frac{{\bf H}^{-1} {\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\bf H}^{-1}}{1+{\mathbf {v}}^{{\sf T}} {\bf H}^{-1} {\mathbf {u}}} \end{eqnarray} (24)to get a rank-one update   \begin{eqnarray} \frac{\partial {\mathbf {E}}_{\setminus 1} }{\partial x_n} = -({\mathbf {E}}{\mathbf {u}}{\mathbf {v}}^{{\sf T}} {\mathbf {E}})_{\setminus 1} \end{eqnarray} (25)for a small perturbation such that $${\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}\ll 1$$. To go from a 2-D solution to a 2.5-D solution, a correction term k2T appears in the system matrices so that the shape matrices S(e) are extended to become S(e) + k2T(e) for spatial wave-number k.   \begin{eqnarray} {\mathbf {T}}^{(e)} = \frac{1}{2|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (26)For the 2.5-D position Jacobian, this change adds a new partial derivative term ∂T(e)/∂xn  \begin{eqnarray} \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} \, \rightarrow \, \frac{\partial ({\mathbf {S}}^{(e)}+k^2{\mathbf {T}}^{(e)})}{\partial x_n} = \frac{\partial {\mathbf {S}}^{(e)}}{\partial x_n} + k^2 \frac{\partial {\mathbf {T}}^{(e)}}{\partial x_n} \end{eqnarray} (27)where ∂S(e)/∂xn is already available from the 2-D calculations, and the additional term ∂T(e)/∂xn may be derived   \begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{\partial |\det {\mathbf {E}}|^{-1}}{\partial x_n} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (28)but we already have $$\partial |\det {\mathbf {E}}|^{-1} / \partial x_n$$ from deriving the partial derivatives of the S(e) term (23) giving   \begin{eqnarray} \frac{\partial {\mathbf {T}}^{(e)} }{\partial x_n} = \frac{1}{2}\frac{ {\mathbf {v}}^{{\sf T}} {\mathbf {E}}{\mathbf {u}}}{|\det {\mathbf {E}}|} \left[\begin{array}{c{@}{\quad }c{@}{\quad }c}2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right] \frac{1}{12} \end{eqnarray} (29)for a linear interpolatory shape function E perturbing a node at row u and column v. This adjoint or rank-one perturbation method for the 2.5-D electrode position Jacobian may be calculated much more quickly than a direct perturbation method because the system matrices do not need to be recalculated and inverted to determine the change in measurements due to electrode movement. The 2.5-D electrode position Jacobian (13) was found to be 25.9 times faster than the equivalent 3-D rank-one update method (Gómez-Laberge & Adler 2008), implementing (12) in 3-D, for mesh geometries used in this work (Intel Core i5-2500K 4-core processor at 3.30 GHz with 32 GB memory). 2-D electrode position Jacobian estimates differ significantly from 3-D solutions (Fig. 8), so have not been presented in Fig. 9. The computational cost of the 2-D rank-one update method for calculating the electrode position Jacobian (Gómez-Laberge & Adler 2008) was orders of magnitude faster than a naïve 2-D perturbation method. The 2-D rank-one update was 7.1 times faster than the 2.5-D method across most mesh sizes, which is accounted for by the numerical integration implied by (13). There are likely to be further gains from optimizing this implementation for multiple processing cores because key portions of the Jacobian calculation (14), (15), (26), (29) can be performed in parallel and the Jacobian typically consumes a significant portion of the total calculation time in each Gauss–Newton iteration (Boyle et al. 2012b). Mesh density was measured as the inverse of the maximum element height h (elements per metre) for both 2-D and 3-D meshes. The relative speed-up for a particular mesh density 1/h grows as a function of the number of FEM mesh nodes n and elements which must be calculated in the Jacobian where n = O(h−2) in 2-D and n = O(h−3) in 3-D; the larger the mesh the greater the benefit conferred by the 2.5-D Jacobian approach. Figure 9. View largeDownload slide The 2.5-D Jacobian calculations scale with FEM node count: the 2.5-D movement Jacobian speed advantage over a rank-one 3-D calculation (Gómez-Laberge & Adler 2008) grows as mesh density (FEM elements per metre) 1/h, where h is the maximum element height for the entire mesh. Error bars show max/min run times over 20 runs, the run time is closely related to the number of nodes in the FEM mesh where 2-D meshes have far fewer nodes for the same mesh density. Figure 9. View largeDownload slide The 2.5-D Jacobian calculations scale with FEM node count: the 2.5-D movement Jacobian speed advantage over a rank-one 3-D calculation (Gómez-Laberge & Adler 2008) grows as mesh density (FEM elements per metre) 1/h, where h is the maximum element height for the entire mesh. Error bars show max/min run times over 20 runs, the run time is closely related to the number of nodes in the FEM mesh where 2-D meshes have far fewer nodes for the same mesh density. 6 RESULTS The column ℓ2-norm sensitivity (the diagonal of $${{\bf J}}^{{\sf T}} {{\bf J}}$$) was plotted by replacing reconstructed resistivity with the log of estimated sensitivity. Sensitivity in these plots was expected to be greatest near the electrodes and diminish elsewhere. Simple Dirichlet boundary conditions away from the electrodes, used in these simulations, introduced errors, which were observable as variations in sensitivity at unexpected locations. We wished to model an approximately half-space forward model but unexpected increases in sensitivity near the sides and bottom were found to be caused by the boundary conditions which were deflecting current flow. Boundary condition errors can be corrected in a number of ways, the simplest of which is to increase the modelled domain until the error is small enough. One could, alternatively, implement appropriate ‘infinite elements’ at the boundary (Babus̆ka 1972). Another method is to estimate a ‘primary’ field for each stimulus using an analytic half-space model or very detailed one-time-use FEM mesh, and then calculate a ‘secondary’ field update on a smaller FEM with different resistivity as a correction (Günther et al. 2006). The primary-secondary type methods rely on small changes in resistivity far enough from the boundary to leave the ‘primary’ field largely unperturbed and it is not immediately obvious how this method may be affected by electrode movements or surface deformation without recalculating the primary field at each update. Neumann and mixed boundary conditions away from the electrodes were not considered. We have used an expanded model, in the interests of reliable results under deformed boundaries, at the expense of some lost computational efficiencies. Regions of high sensitivity were initially noted at depth where little sensitivity was expected. To determine how far the FEM model boundaries needed to be extended, an analytic PEM half-space model was compared to CEM homogeneous resistivity FEM simulated measurements. The model boundaries were extended approximately one electrode array length in each of the +x, −x and −z directions. The boundary extension reduced boundary condition related errors in simulated measurements to within measured noise levels and removed the artefacts from the sensitivity plots. The resistivity sensitivity S was plotted relative to the maximum sensitivity,   \begin{eqnarray} {\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V} \end{eqnarray} (30)using the conductivity Jacobian Jσ, measurement inverse covariance/weighting W, and element volumes $$\bf V$$, for line#1 (March 2008) and line#5 (February 2013) using the surveyed locations (Fig. 10). The region near the electrodes has been presented with annotations matching Fig. 6, as well as an image of the complete model. The initial and final resistivities were independently reconstructed using the surveyed locations (Fig. 11) and the difference between initial and final was used to create the expected resistivity change (Figs 12e and f). The initial resistivity for line#1 closely matched those published in Wilkinson et al. (2010) (Fig. 6) and achieved a similar <1 per cent RMS measurement misfit relative to a homogeneous resistivity model. Qualitatively and quantitatively, the two reconstructions (Figs 6 and 11a) are very similar. Figure 10. View largeDownload slide Inverse model parametrization with colours showing relative sensitivity $${\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V}$$ as S/max (S) for homogeneous resistivity; (a) line#1 March 2008, and (b) line#5 February 2013; note the distinct difference in slope profile between (a) and (b). Figure 10. View largeDownload slide Inverse model parametrization with colours showing relative sensitivity $${\bf S} = {\rm diag}({{\bf J}}_\sigma ^{{\sf T}} {\mathbf {W}}{{\bf J}}_\sigma )/{\bf V}$$ as S/max (S) for homogeneous resistivity; (a) line#1 March 2008, and (b) line#5 February 2013; note the distinct difference in slope profile between (a) and (b). Figure 11. View largeDownload slide Resistivity-only reconstructions using true electrode positions as measured by RTK GPS; (a) line#1 March 2008 to (b) April 2009 (λ = 52.1, σ0 = 31.3 Ωm), and (c) line#5 February 2013 to (d) February 2014 (λ = 42.7, σ0 = 25.6 Ωm). Figure 11. View largeDownload slide Resistivity-only reconstructions using true electrode positions as measured by RTK GPS; (a) line#1 March 2008 to (b) April 2009 (λ = 52.1, σ0 = 31.3 Ωm), and (c) line#5 February 2013 to (d) February 2014 (λ = 42.7, σ0 = 25.6 Ωm). Figure 12. View largeDownload slide Change in resistivity and electrode movement for joint movement reconstructions (λσ = 0.1, λx = 0.07), (a–d) line#1, March 2008 to April 2009, and (e–h) line#5, February 2013 to February 2014. Figure 12. View largeDownload slide Change in resistivity and electrode movement for joint movement reconstructions (λσ = 0.1, λx = 0.07), (a–d) line#1, March 2008 to April 2009, and (e–h) line#5, February 2013 to February 2014. Electrode movements were initially reconstructed without allowing resistivity change using an independent implementation of the electrode movement iterative solver. The observed behaviour of the algorithm (Fig. 13) was to first minimize the error in the electrode spacing (corresponding to the largest measurement misfit), and then to approach a more ‘correct’ solution. Artefacts of this approach to the minimum over the optimization surface remain where large steps in the true electrode movement are reconstructed as a balanced step by adjacent electrodes that come close to the true displacement between electrodes near that movement. The joint inversion code was checked against this result by setting the movement-resistivity balance parameter β to strongly favour electrode movement. Reconstructions showed no resistivity change and movements that were in close agreement with the movement-only reconstruction code. Small variations still existed between the two results due to differences in the Gauss–Newton implementation and inexact line search. These variations were small with respect to the overall electrode movement solution. When electrode movements were reconstructed with resistivity changes (Fig. 12), some portion of the reconstructed electrode movement was lost in favour of reconstructed resistivity change. Figure 13. View largeDownload slide Electrode movement without allowing for resistivity changes, iterations for (a) line#1, March 2008 to April 2009, and (b) line#5 movement, February 2013 to February 2014. Figure 13. View largeDownload slide Electrode movement without allowing for resistivity changes, iterations for (a) line#1, March 2008 to April 2009, and (b) line#5 movement, February 2013 to February 2014. Due to the large electrode movements, it was found helpful to perform a crude version of successive relaxation. The first three iterations of the Gauss–Newton reconstruction were performed with an electrode movement hyperparameter that was an order of magnitude larger than following iterations. Without this adaptation, the reconstructed electrode movements showed poor agreement with measured locations, presumably because the Gauss–Newton iterations were trapped in a local minimum which favoured constructing resistivity artefacts near the electrodes. Exploring the space of hyperparameters near the selected hyperparameter did not reveal one which achieved better electrode movement reconstruction. 7 DISCUSSION Resistivity was reconstructed for measured initial and final electrode locations (Fig. 11) which serve as an ‘ideal’ reconstruction. Resistivity and electrode displacement were simultaneously reconstructed for a survey located on a slowly moving landslide. Results exhibit some measure of oscillatory artefacts in the reconstructed movement (Figs 12c and g). The resistivity distribution for line#1 (March 2008–April 2009; Figs 11a and b) changed by a relatively small amount when using true electrode locations before and after movement. This would suggest that, beyond the ground motion at the surface, no structural changes in the near surface seem to have occurred. It seems plausible that the increased area of low resistivity WMF might be indicative of increased saturation of the soil which led to the translational slide of WMF material moving over SSF substrate at the surface. The resistivity changes for line#5 (February 2013–February 2014; Figs 11c and d), show a significantly different distribution after ground movement, which is interesting given that the two electrode lines are within 40 m of each other. The line#5 measurements (Fig. 11d) occurred after a very wet summer and winter period where there was a lot of seepage at the base of the lobe causing the deeper reductions in resistivity between z = 40 and 60 m. At the surface of the lobes, resistivity increased due to cracking of the top layer. A cracked surface experienced accelerated evaporation over increased surface area, resulting in localized resistivity increase. To a lesser extent, areas affected by surface cracking also showed increased resistivity due to the change in topography: relatively little current would be conducted across the air gap in cracks, contributing to an average increase in bulk conductivity, but the cracks do increase the surface along which current flows resulting in an effective increase in resistivity. We take particular note of the change in SSF resistivity around x = 80 m which may have developed vertical connectivity between the overlying WMF and RMF below, allowing vertical drainage. The flow would be from the saturated WMF, along the WMF-SSF boundary to the surface, then to a region of vertical connectivity downward through the SSF (x = 80 m), and then into the RMF where it has pooled underground. This proposed flow path might also explain the increase in resistivity in the SSF at x = 60 m: if the vertical connectivity were in a roughly vertical plane, it would cut off the outer section of the SSF and that outer section would drain downwards into the RMF leading to an increase in resistivity. The deeper segment of the SSF (x > 80 m) would maintain its resistivity because the general connectivity and saturation have not changed by much. The resistivity change may also be induced by model error in electrode placement or topography: the 2.5-D model limits fidelity in some respects. The poor quality of the line#5 post-movement data, as measured by reciprocal error, may be the cause of these changes, though the locations of the reciprocal errors along the electrode array were distributed along the length of the array so that we expect no concentrated region of low sensitivity that could cause resistivity changes in the reconstruction such as those observed Fig. 11(d). Given this speculation, it would be interesting to investigate this potential vertical fault in the SSF. Perhaps it is an indication of a major ground movement, still to come, as the lower slope drainage has changed significantly. The change in drainage may also help to stabilize the lower slope by providing a drainage path at-depth which will allow surface material in the SSF to consolidate. This stabilizing effect has been observed on other nearby slopes in previous years. It has been postulated in Uhlemann et al. (2017) that reactivation of the slope at line#1 was stabilized due to slope movement which caused preferential flow paths to open, lowering pore pressures on the slip surface of the lobe, thereby stabilizing the lobe. For a half-space model with a homogeneous resistivity, electrode positions are not unique. A translation of the entire set of electrodes will give identical measurements. Likewise, a scaling of all electrode positions is equivalent to a change in the homogeneous resistivity. When conductivities are inhomogeneous the electrode locations are somewhat fixed by the locations of the inhomogeneities. Examples of electrode position non-uniqueness manifested itself in this data as large oscillations in the reconstructed electrode movement when no measures were taken to address the issue. Fixing the location of three electrodes at the upslope and downslope ends of the electrode array (Figs 12c and g) nearly eliminated these oscillations. As seen in the line#5 data, this is not necessarily a correct assumption, as both the top and bottom of a landslide may move, leading to resistivity artefacts. We infer that fixing these electrode locations was sufficient to damp the reconstructed movement’s oscillations because it fixes the relationship between a stimulus current, a measured voltage, and electrode separation (distance). Smoothing-type regularization of resistivity, used in this work, then controls how strictly the selected scaling of electrode separation is enforced. For example, electrodes that have contracted together in a region might cause a conductive artefact to be reconstructed under those electrodes. Increasing the resistivity regularization may suppress these artefacts and cause movement to be reconstructed by contracting the local electrode spacing to account for smaller than expected voltage measurements in the region. Both unscaled and log scaled electrode movement were tested and found to give results with similar absolute positional error. We have elected to present the unscaled electrode movement in our reconstructions as it is a less restrictive choice. One could imagine low angle slopes, wetlands or floodplains where the expected direction of movement would not be known a priori. Peat wetlands, muskeg, and most ground exposed to deep frost experience seasonal expansion and contraction, due to freeze-thaw cycles, water table changes, or water and gas accumulation and evaporation, which can result in uplift and ground shifting in directions other than downslope (Taber 1930; Hansell et al. 1983; Price 2003; Strack et al. 2006; Uhlemann et al. 2016c). A further reason to avoid dependence on the log movement scaling is the extension of this work to transverse movements where the restriction to movements only to one side of the array seems inappropriate. In the data sets examined here, there are transverse movements which were caused by material accumulating at the toe of the landslide and towards the edges of the earth flow. These transverse movements can be predicted for this particular data set based on the pre-existing topology: line#1 electrodes moved east, downhill into a gully, while the line#5 electrodes moved west, downhill into the same gully. Reconstructions for line#1 generally matched the true electrode locations within 0.20 m for movements of up to 1.46 m with the exception of electrode #9 and the three electrodes #6–#8 at the step in electrode position (Fig. 12a). Compared to Wilkinson et al. (2010) (≤0.2 m position error), these results are marginally less accurate. It seems probable that our results for the line#1 data differ from those of Wilkinson et al. (2010) due to the simultaneous resistivity and electrode position reconstruction, presented here, which removes artificial ordering constraints that occur with the sequential method of Wilkinson et al. (2010). Another source of differences in our results with respect to Wilkinson et al. (2010) is the restriction to downslope movements by a log parametrization in Wilkinson et al. (2010). As mentioned previously, trials of this log parametrization method in our simultaneous resistivity and electrode position inversion did not lead to improved results. Our results for the line#5 data appear to closely correlate with Wilkinson et al. (2016), where reconstructed electrode locations were generally within 0.2 m excepting some electrodes with errors up to 1.0 m, though results are not presented in as much detail in that case. In contrast to Wilkinson et al. (2016), our reconstructions do not address movements transverse to the electrode line. Movement reconstructions for line #5 do not appear to be particularly accurate, perhaps due to the more significant resistivity changes inferred in the reconstruction and more widespread translational failure of the slope which shifted electrodes over most of the resistivity structure (Fig. 12g). These might be addressed by identifying the covariance between movement and resistivity change within the joint reconstruction algorithm regularization. It is also possible that some of the error in reconstructed electrode position may be due to the FEM discretization. An approach such as the Fréchet method for electrode movement may help to alleviate such issues, though in general it produces the same solutions as a 3-D rank-one update method (Dardé et al. 2012; Boyle et al. 2017). Adjusting the relationship between resistivity and movement regularization β caused greater electrode displacement error as resistivity regularization was reduced. These movement magnitudes represented movement of up to 32 per cent of the average 4.73 m electrode spacing, exceeding the joint resistivity-movement methods of Soleimani et al. (2006) which was limited to movements of approximately 1 per cent of electrode spacing. 8 CONCLUSIONS This work demonstrated the practical application of a joint electrode movement-resistivity reconstruction using an iterative Gauss–Newton regularized solver. The electrode position Jacobian was calculated on the current resistivity at each iteration. Reconstructions show reasonable agreement with RTK GPS measured electrode locations, available resistivity estimates and geological structure. The initial reconstructed resistivity model, used as a starting point for the electrode movement and resistivity reconstruction, was in close agreement with prior work (Figs 6 and 11a) (Wilkinson et al. 2010). Reconstructed changes in resistivity (Fig. 12) showed considerable variation, particularly around the region at the toe of the landslide. These changes in resistivity could be indicative of water saturation changes due to water seepage at the toe of the landslide or other geological causes. Another possibility is that the resistivity changes represent artefacts due to transverse and normal components of electrode movement which were not accounted for in these reconstructions. We note that, in general, even when electrode displacements were not entirely accurate compared to true electrode positions, the error in the estimated and true distances between electrode positions was quite accurate after one or two Gauss–Newton iterations. Errors in electrode spacing were distributed fairly evenly across the electrode array, after which the displacements shifted towards their true positions in most cases. This suggests that a parametrization for electrode movement that encompasses electrode spacing may lead to improved outcomes. The results suggest that electrode grids are effective not only for resistivity monitoring, but also as a means of ground motion detection which may provide a cost-effective approach for landslide monitoring. We are encouraged by these results and expect that with new protocols which measure between electrode lines, higher quality electrode position reconstructions will be possible which incorporate normal, lateral, and transverse electrode movements, as observed in the data sets presented here. Acknowledgements Identification of the Hollin Hill site, geology, field work and electrical measurements were carried out by the British Geological Survey. Thanks go out to the researchers at the British Geological Survey, Geophysical Tomography Team for generously sharing their time and data. They were supported by the Natural Environment Research Council (NERC) and their contributions to this work are published with the permission of the Executive Director of the British Geological Survey. In fond memory of Steve Gibson, and with thanks to Josie Gibson. Code for the algorithms presented in this work extend EIDORS, an open source Matlab toolkit for impedance imaging inverse problems (Adler & Lionheart 2006; Adler et al. 2015, 2017), and use NetGen for meshing (Schöberl 1997). This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC). Footnotes 1 Leica System 1200 RTK GPS in ‘kinematic mode’ real-time correction achieves as much as 10 mm (Root Mean Squared (RMS)) horizontal and 20 mm vertical accuracies (Merritt 2014). 2 Longitudinal movement being movement inline with the electrodes and along the surface. 3 Matlab interp1(X,Y,Xq,‘pchip’) 1-D interpolant based on Hermite derivatives. REFERENCES Adler A., Lionheart W.R.B., 2006. Uses and abuses of EIDORS: an extensible software base for EIT, Physiol. Meas. , 27( 5), S25– S42. https://doi.org/10.1088/0967-3334/27/5/S03 Google Scholar CrossRef Search ADS PubMed  Adler A., Guardo R., Berthiaume Y., 1996. Impedance imaging of lung ventilation: do we need to account for chest expansion?, IEEE Trans. Biomed. Eng. , 43( 4), 414– 420. https://doi.org/10.1109/10.486261 Google Scholar CrossRef Search ADS PubMed  Adler A., Gaburro R., Lionheart W. R.B., 2011. Electrical impedance tomography, in Mathematical Methods in Imaging , chap. 14, pp. 601– 654, ed. Scherzer O., Springer Science+Business Media. Adler A., Boyle A., Crabb M.G., Gagnon H., Grychtol B., Lesparre N., Lionheart W.R.B., 2015. EIDORS Version 3.8, in International Conference on Biomedical Applications of Electrical Impedance Tomography , Neuchâtel, Switzerland. Adler A., Boyle A., Braun F., Crabb M.G., Grychtol B., Lionheart W. R.B., Tregidgo H.F.J., Yerworth R., 2017. EIDORS 3.9, in 18th International Conference on Biomedical Applications of Electrical Impedance Tomography , Dartmouth, USA. Babus̆ka I., 1972. The finite element method for infinite domains. I, Math. Comput. , 26( 117), 1– 11. https://doi.org/10.1090/S0025-5718-1972-0298969-2 Google Scholar CrossRef Search ADS   Boyle A., 2010. The effect of boundary shape deformation on two-dimensional electrical impedance tomography, Master’s thesis , Carleton University, Ottawa, Canada. Boyle A., 2016. Geophysical applications of electrical impedance tomography, PhD thesis , Carleton University, Ottawa, Canada. Boyle A., Adler A., 2010. Electrode models under shape deformation in electrical impedance tomography, in 11th Conf. Electrical Impedance Tomography , University of Florida, Gainesville, USA. Google Scholar CrossRef Search ADS   Boyle A., Adler A., 2011. Impact of electrode area, contact impedance and boundary shape on EIT images, Physiol. Meas. , 32( 7), 745– 754. https://doi.org/10.1088/0967-3334/32/7/S02 Google Scholar CrossRef Search ADS PubMed  Boyle A., Adler A., 2016. Modelling with 2.5D approximations, in 17th Conference on Electrical Impedance Tomography , Stockholm, Sweden. Boyle A., Adler A., Lionheart W.R.B., 2012a. Shape deformation in two-dimensional electrical impedance tomography, IEEE Trans. Med. Imaging , 31( 12), 2185– 2193. https://doi.org/10.1109/TMI.2012.2204438 Google Scholar CrossRef Search ADS   Boyle A., Borsic A., Adler A., 2012b. Addressing the computational cost of large EIT solutions, Physiol. Meas. , 33( 5), 787– 800. https://doi.org/10.1088/0967-3334/33/5/787 Google Scholar CrossRef Search ADS   Boyle A., Wilkinson P., Chambers J., Lesparre N., Adler A., 2014. Slope stability monitoring through impedance imaging, in 15th Conf. Electrical Impedance Tomography , Gananoque, Canada. Boyle A., Crabb M.G., Jehl M., Lionheart W.R.B., Adler A., 2017. Methods for calculating the electrode movement Jacobian for impedance imaging, Physiol. Meas. , 38( 3), 555– 574. https://doi.org/10.1088/1361-6579/aa5b78 Google Scholar CrossRef Search ADS PubMed  Carlson R., Fritsch F., 1989. An algorithm for monotone piecewise bicubic interpolation, SIAM J. Numer. Anal. , 26( 1), 230– 238. https://doi.org/10.1137/0726013 Google Scholar CrossRef Search ADS   Chambers J. et al.  , 2011. Three-dimensional geophysical anatomy of an active landslide in Lias Group mudrocks, Cleveland Basin, UK, Geomorphology , 125( 4), 472– 484. https://doi.org/10.1016/j.geomorph.2010.09.017 Google Scholar CrossRef Search ADS   Cormen T., Leiserson C., Rivest R., Stein C., 1990. Introduction to Algorithms , 3rd edn, The MIT Press. Dardé J., Hakula H., Hyvönen N., Staboulis S., 2012. Fine-tuning electrode information in electrical impedance tomography, Inverse Probl. Imaging , 6( 3), 399– 421. https://doi.org/10.3934/ipi.2012.6.399 Google Scholar CrossRef Search ADS   Dardé J., Hyvönen N., Seppänen A., Staboulis S., 2013a. Simultaneous reconstruction of outer boundary shape and admittivity distribution in electrical impedance tomography, SIAM J. Imaging Sci. , 6( 1), 176– 198. https://doi.org/10.1137/120877301 Google Scholar CrossRef Search ADS   Dardé J., Hyvönen N., Seppänen A., Staboulis S., 2013b. Simultaneous recovery of admittivity and body shape in electrical impedance tomography: an experimental evaluation, Inverse Probl. , 29( 8), 1– 16. Google Scholar CrossRef Search ADS   Dey A., Morrison H., 1979. Resistivity modelling for arbitrarily shaped two-dimensional structures, Geophys. Prospect. , 27, 106– 136. https://doi.org/10.1111/j.1365-2478.1979.tb00961.x Google Scholar CrossRef Search ADS   Friedel S., Thielen A., Springman S., 2006. Investigation of a slope endangered by rainfall-induced landslides using 3D resistivity tomography and geotechnical testing, J. appl. Geophys. , 60( 2), 100– 114. https://doi.org/10.1016/j.jappgeo.2006.01.001 Google Scholar CrossRef Search ADS   Gómez-Laberge C., Adler A., 2008. Direct EIT Jacobian calculations for conductivity change and electrode movement, Physiol. Meas. , 29( 6), S89– S99. https://doi.org/10.1088/0967-3334/29/6/S08 Google Scholar CrossRef Search ADS PubMed  Günther T., Rücker C., Spitzer K., 2006. Three-dimensional modelling and inversion of DC resistivity data incorporating topography - II. Inversion, Geophys. J. Int. , 166( 2), 506– 517. https://doi.org/10.1111/j.1365-246X.2006.03011.x Google Scholar CrossRef Search ADS   Hansell R., Scott P., Staniforth R., Svoboda J., 1983. Permafrost development in the intertidal zone at Churchill, Manitoba: a possible mechanism for accelerated beach uplift, Arctic , 36( 2), 198– 203. https://doi.org/10.14430/arctic2263 Google Scholar CrossRef Search ADS   Hyvönen N., Seppänen A., Staboulis S., 2014. Optimizing electrode positions in electrical impedance tomography, SIAM J. Appl. Math. , 74( 6), 1831– 1851. https://doi.org/10.1137/140966174 Google Scholar CrossRef Search ADS   Jehl M., Avery J., Malone E., Holder D., Betcke T., 2015. Correcting electrode modelling errors in EIT on realistic 3D head models, Physiol. Meas. , 36( 12), 2423– 2442. https://doi.org/10.1088/0967-3334/36/12/2423 Google Scholar CrossRef Search ADS PubMed  Jomard H., Lebourg T., Binet S., Tric E., Hernandez M., 2007. Characterization of an internal slope movement structure by hydrogeophysical surveying, Terra Nova , 19( 1), 48– 57. https://doi.org/10.1111/j.1365-3121.2006.00712.x Google Scholar CrossRef Search ADS   Jongmans D., Garambois S., 2007. Geophysical investigation of landslides: a review, Bull. Soc. France , 178( 2), 101– 112. https://doi.org/10.2113/gssgfbull.178.2.101 Google Scholar CrossRef Search ADS   Kim J., Yi M., Supper R., Ottowitz D., 2014. Simultaneous inversion of resistivity structure and electrode locations in ERT, in 20th European Meeting of Environmental and Engineering Geophysics: Near Surface 2014 , pp. 560– 564, Athens, Greece. Kuras O., Pritchard J., Meldrum P., Chambers J., Wilkinson P., Ogilvy R., Wealthall G., 2009. Monitoring hyrdaulic processes with automated time-lapse electrical resistivity tomography (ALERT), C. R. Geosci. , 341( 10–11), 868– 885. https://doi.org/10.1016/j.crte.2009.07.010 Google Scholar CrossRef Search ADS   Lapenna V., Lorenzo P., Perrone A., Piscitelli S., Rizzo E., Sdao F., 2005. 2D electrical resistivity imaging of some complex landslides in Lucanian Apennine chain, southern Italy, Geophysics , 70( 3), B11– B18. https://doi.org/10.1190/1.1926571 Google Scholar CrossRef Search ADS   Lebourg T., Binet S., Tric E., Jomard H., Bedoui S.E., 2005. Geophysical survey to estimate the 3D sliding surface and the 4D evolution of the water pressure on part of a deep seated landslide, Terra Nova , 17( 5), 399– 406. https://doi.org/10.1111/j.1365-3121.2005.00623.x Google Scholar CrossRef Search ADS   Lebourg T., Hernandez M., Zerathe S., Bedoui S.E., Jomard H., Fresia B., 2010. Landslides triggered factors analysed by time lapse electrical survey and multidimensional statistical approach, Eng. Geol. , 114( 3–4), 238– 250. https://doi.org/10.1016/j.enggeo.2010.05.001 Google Scholar CrossRef Search ADS   Loke M., 2017. RES2DINV: Rapid 2-D Resistivity & IP Inversion using the Least-squares Method , Geotomosoft Solutions. Loke M., Wilkinson P., Chambers J., 2015. Rapid inversion of data from 2-D and from 3-D resistivity surveys with shifted electrodes, in Near Surface Geoscience 2015-21st European Meeting of Environmental and Engineering Geophysics , Turin, Italy. Loke M., Wilkinson P., Chambers J., 2016. 3-D resistivity inversion with electrodes displacements, in 25th International Geophysical Conference and Exhibition , pp. 806– 810, Adelaide, Australia. Google Scholar CrossRef Search ADS   Loke M., Wilkinson P., Chambers J., Meldrum P., 2017. Rapid inversion of data from 2D resistivity surveys with electrode displacements, Geophys. Prospect. , doi:10.1111/1365-2478.12522. Merritt A., 2014. 4D geophysical monitoring of hydrogeological precursors to landslide activation, PhD thesis , University of Leeds, Leeds, United Kingdom. Merritt A. et al.  , 2014. 3D ground model development for an active landslide in lias mudrocks using geophysical, remote sensing and geotechnical methods, Landslides , 11( 4), 537– 550. https://doi.org/10.1007/s10346-013-0409-1 Google Scholar CrossRef Search ADS   Naudetb V., Lazzari M., Perrone A., Loperte A., Piscitelli S., Lapenna V., 2008. Integrated geophysical and geomorphological approach to investigate the snowmelt-triggered landslide of Bosco Piccolo village (Basilicata, southern Italy), Eng. Geol. , 98( 3–4), 156– 167. https://doi.org/10.1016/j.enggeo.2008.02.008 Google Scholar CrossRef Search ADS   Nocedal J., Wright S., 1999. Numerical Optimization , Springer Science Business Media. Google Scholar CrossRef Search ADS   Oldenborger G., Routh P., Knoll M., 2005. Sensitivity of electrical resistivity tomography data to electrode position errors, Geophys. J. Int. , 163( 1), 1– 9. https://doi.org/10.1111/j.1365-246X.2005.02714.x Google Scholar CrossRef Search ADS   Perrone A., Iannuzzi A., Lapenna V., Lorenzo P., Piscitelli S., Rizzo E., Sdao F., 2004. High-resolution electrical imaging of the Varco d’Izzo earthflow (southern Italy), J. Appl. Geophys. , 56( 1), 17– 29. https://doi.org/10.1016/j.jappgeo.2004.03.004 Google Scholar CrossRef Search ADS   Piegaria E., Cataudella V., Maio R.D., Milano L., Nicodemi M., Soldovieri M., 2009. Electrical resistivity tomography and statistical analysis in landslide modelling: a conceptual approach, J. Appl. Geophys. , 68( 2), 151– 158. https://doi.org/10.1016/j.jappgeo.2008.10.014 Google Scholar CrossRef Search ADS   Price J., 2003. Role and character of seasonal peat soil deformation on the hydrology of undisturbed and cutover peatlands, Water Resour. Res. , 39( 9), SBH–1– SBH– 10. Rücker C., Günther T., 2011. The simulation of finite ERT electrodes using the complete electrode model, Geophysics , 76( 4), F227– F238. https://doi.org/10.1190/1.3581356 Google Scholar CrossRef Search ADS   Sass O., Bell R., Glade T., 2008. Comparison of GPR, 2D-resistivity and traditional techniques for the subsurface exploration of the Öschingen landslide, Swabian Alb (Germany), Geomorphology , 93( 1-2), 89– 103. https://doi.org/10.1016/j.geomorph.2006.12.019 Google Scholar CrossRef Search ADS   Schöberl J., 1997. NETGEN an advancing front 2D/3D-mesh generator based on abstract rules, Comput. Vis. Sci. , 1( 1), 41– 52. https://doi.org/10.1007/s007910050004 Google Scholar CrossRef Search ADS   Soleimani M., Gómez-Laberge C., Adler A., 2006. Imaging of conductivity changes and electrode movement in EIT, Physiol. Meas. , 27( 5), S103– S113. https://doi.org/10.1088/0967-3334/27/5/S09 Google Scholar CrossRef Search ADS PubMed  Somersalo E., Cheney M., Isaacson D., 1992. Existence and uniqueness for electrode models for electric current computed tomography, SIAM J. Appl. Math. , 52( 4), 1023– 1040. https://doi.org/10.1137/0152060 Google Scholar CrossRef Search ADS   Strack M., Kellner E., Waddington J., 2006. Effect of entrapped gas on peatland surface level fluctuations, Hydrol. Process. , 20( 17), 3611– 3622. https://doi.org/10.1002/hyp.6518 Google Scholar CrossRef Search ADS   Supper R. et al.  , 2014. Geoelectrical monitoring: an innovative method to supplement landslide surveillance and early warning, Near Surf. Geophys. , 12( 1), 133– 150. Suzuki K., Higashi S., 2001. Groundwater flow after heavy rain in landslide-slope area from 2-D inversion of resistivity monitoring data, Geophysics , 66( 3), 733– 743. https://doi.org/10.1190/1.1444963 Google Scholar CrossRef Search ADS   Taber S., 1930. The mechanics of frost heaving, J. Geol. , 38( 4), 303– 317. https://doi.org/10.1086/623720 Google Scholar CrossRef Search ADS   Uhlemann S., Hagedorn S., Dashwood B., Maurer H., Gunn D., Dijkstra T., Chambers J., 2016a. Landslide characterization using P- and S-wave seismic refraction tomography — the importance of elastic moduli, J. Appl. Geophys. , 134( 1), 64– 76. https://doi.org/10.1016/j.jappgeo.2016.08.014 Google Scholar CrossRef Search ADS   Uhlemann S. et al.  , 2016b. Assessment of ground-based monitoring techniques applied to landslide investigations, Geomorphology , 253( 1), 438– 451. https://doi.org/10.1016/j.geomorph.2015.10.027 Google Scholar CrossRef Search ADS   Uhlemann S., Sorensen J., House A., Wilkinson P., Roberts C., Gooddy D., Binley A., Chambers J., 2016c. Integrated time-lapse geoelectrical imaging of wetland hydrological processes, Water Resour. Res. , 52( 3), 1607– 1625. https://doi.org/10.1002/2015WR017932 Google Scholar CrossRef Search ADS   Uhlemann S. et al.  , 2017. 4D imaging of moisture dynamics during landslide reactivation, J. geophys. Res. , 122( 1), 398– 418. https://doi.org/10.1002/2016JF003983 Google Scholar CrossRef Search ADS   Wagner F.M., Bergmann P., Rücker C., Wiese B., Labitzke T., Schmidt-Hattenberger C., Maurer H., 2015. Impact and mitigation of borehole related effects in permanent crosshole resistivity imaging: an example from the Ketzin CO2 storage site, J. Appl. Geophys. , 123( 1), 102– 111. https://doi.org/10.1016/j.jappgeo.2015.10.005 Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Lelliot M., Wealthall G., Ogilvy R., 2008. Extreme sensitivity of crosshole electrical resistivity tomography measurements to geometric errors, Geophys. J. Int. , 173( 1), 49– 62. https://doi.org/10.1111/j.1365-246X.2008.03725.x Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Meldrum P., Gunn D., Ogilvy R., Kuras O., 2010. Predicting the movements of permanently installed electrodes on an active landslide using time-lapse geoelectrical resistivity data only, Geophys. J. Int. , 183( 2), 543– 556. https://doi.org/10.1111/j.1365-246X.2010.04760.x Google Scholar CrossRef Search ADS   Wilkinson P., Chambers J., Uhlemann S., Meldrum P., Smith A., Dixon N., Loke M., 2016. Reconstruction of landslide movements by inversion of 4-D electrical resistivity tomography monitoring data, Geophys. Res. Lett. , 43( 3), 1166– 1174. https://doi.org/10.1002/2015GL067494 Google Scholar CrossRef Search ADS   Zhou B., Dhalin T., 2003. Properties and effects of measurement errors on 2D resistivity imaging surveying, Near Surf. Geophys. , 1( 3), 105– 117. https://doi.org/10.3997/1873-0604.2003001 © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 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. ### DeepDyve Freelancer ### DeepDyve Pro Price FREE$49/month

\$360/year
Save searches from Google Scholar, PubMed
Create lists to organize your research
Export lists, citations