A linked geomorphological and geophysical modelling methodology applied to an active landslide

Moisture-induced landslides are a global geohazard; mitigating the risk posed by landslides requires an understanding of the hydrological and geological conditions present within a given slope. Recently, numerous geophysical studies have been attempted to characterise slow-moving landslides, with an emphasis on developing geoelectrical methods as a hydrological monitoring tool. However, landslides pose specific challenges for processing geoelectrical data in long-term monitoring contexts as the sensor arrays can move with slope movements. Here we present an approach for processing long-term (over 8 years) geoelectrical monitoring data from an active slow-moving landslide, Hollin Hill, situated in Lias rocks in the southern Howardian Hills, UK. These slope movements distorted the initial setup of the monitoring array and need to be incorporated into a time-lapse resistivity processing workflow to avoid imaging artefacts. We retrospectively sourced seven digital terrain models to inform the topography of our imaging volumes, which were acquired by either Unmanned Aerial Vehicle (UAV)-based photogrammetry or terrestrial laser ranging systems. An irregular grid of wooden pegs was periodically surveyed with a global position system, from which distortions to the terrain model and electrode positions can be modelled with thin plate splines. In order to effectively model the time-series electrical resistivity images, a baseline constraint is applied within the inversion scheme; the result of the study is a time-lapse series of resistivity volumes which also incorporate slope movements. The workflow presented here should be adaptable for other studies focussed on geophysical/geotechnical monitoring of unstable slopes.


Introduction
Landslides are a global phenomenon, resulting in severe economic and societal losses, and as such represent a significant geohazard. The majority of land slip events are moisture-induced (Gasmo et al. 2000), whereby increases in subsurface moisture change the pore pressure conditions which consequently affect the shear strength within a slope, resulting in slope failure (e.g., Terzaghi 1936). In order to manage this hazard, it is necessary to characterise landslide bodies both internally and externally. The external geomorphology of unstable slopes can be characterised directly with observations, aerial imagery and laser ranging methods. Determining the internal structure of landslides remains more challenging, often practitioners need to rely on point sensors and physical samples (recovered from pits or core). Over the past few decades, several studies and reviews have investigated the use of geophysical methods for landslide investigation since they are spatially sensitive, non-invasive and comparatively inexpensive relative to conventional shallow borehole investigations (Jongmans and Garambois 2007;Pazzi et al. 2019;Whiteley et al. 2019). Bogoslovsky and Ogilvy (1977) first demonstrated that geoelectrical techniques could be used to make interpretations on the structure of landslides and likely hydrological conditions, as relationships between electrical resistivity and moisture content have been long established (e.g., Archie 1947).
Numerous studies have shown that electrical resistivity imaging (ERI), also known as electrical resistivity tomography (ERT), can be an invaluable aid in interpreting changes in near surface hydrologic conditions (Binley et al. 2015;Brunet et al. 2010;Chambers et al. 2014;Johnson et al. 2017;McLachlan et al. 2020;Perrone et al. 2014;Revil et al. 2020;Uhlemann et al. 2017;Uhlemann et al. 2016b). In the absence of any changes to geological structure, changes in electrical resistivity should be due to changes in temperature and the pore fluid (saturation/salinity) in the subsurface (Waxman and Smits 1968). Hence, ERI has proven to be a powerful tool when used in a hydrological monitoring context (Johnson et al. 2017;Uhlemann et al. 2017). The motivation for conducting time-lapse geoelectrical surveys on landslides is clear; the relationships between moisture content and resistivity show that these methods can be used to infer the hydrological state of a hillslope and by extension shear strength and liquid limits, key parameters in estimating slope stability. For this reason, the number of geoelectrical studies in landslide prone areas has been increasing in recent years (Pazzi et al. 2019;Whiteley et al. 2019). Uhlemann et al. (2017) investigated the use of the Waxman-Smits relationship (Waxman and Smits 1968) for monitoring seasonal moisture content fluctuations in an active landslide over a 3-year time period, showing that elevated moisture content derived from ERI measurements can be associated with slope movements. Crawford and Bryson (2018) presented a novel study directly relating electrical conductivity (the inverse of resistivity) to soil suction which is then used to compute an unsaturated shear strength. Recently, Revil et al. (2020) demonstrated the use of the time domain induced polarisation (IP) method for use in clay rich, landslide prone, materials, transforming their geoelectrical models into both soil moisture content estimates and cation exchange capacities through petrophysical calibration. Once these parameters have been estimated, a volumetric approximation of permeability can be attempted (Soueid , hence this method may have future implications for coupled geoelectrical and hydrological modelling. The focus here is the more widely used ERI method for monitoring landslides. The question addressed in this study is how to process longterm data on an active landslide? It is imperative that workflows are developed to process time-lapse geoelectrical datasets in a timely and robust manner, as by their nature landslides can have multiple data processing challenges associated with them. Misplaced electrodes (in the geophysical model) have potential fields which are incorrectly reproduced in geoelectrical imaging, therefore the user of geoelectrical monitoring must have a good understanding of both the surface topography and electrode placement of a given field site before attempting any geophysical method. Foremost, if the landslide is active, then it is likely that the surface will be altered throughout the monitoring period, and secondly permanently installed electrodes are likely to have translated with surface movements. The latter has been addressed in the literature as the movements of electrodes mask any changes in resistivity due to moisture contents and can cause significant artefacts in the resistivity images if not accounted for in geophysical processing (Uhlemann et al. 2017;Uhlemann et al. 2015;Wilkinson et al. 2016;Wilkinson et al. 2010). Uhlemann et al. (2015) demonstrate three interpolation techniques for interpolating electrode movements from sparse topographic information at a single point in time and Wilkinson et al. (2016) reconstruct landslide movements from 4D ERI monitoring data, whereby changes in the measured transfer resistances are modelled in terms of the electrode displacements. However, neither of the aforementioned methods addresses changes in topography, although Loke et al. (2018) reconstruct topography changes from modelling electrode displacements for 2D monitoring setups. Currently, it is logistically difficult, and cost inefficient, to acquire digital elevation models (DEMs) at a temporal resolution needed for effective geoelectrical monitoring (every 2-3 days in this study). However sparse monitoring of discrete topographic points is more accessible. For example, Le Breton et al. (2019) demonstrate a relatively low-cost monitoring system, where unwrapped phase changes recorded between a radio transmitter and a network of receivers (placed on the moving slope) are translated into one-dimensional movements. The reconstruction of electrode movements with geoelectrical data (Wilkinson et al. 2016) could also be used for this purpose. Here we manually survey gridded markers on a slope with repeated field visits.

Motivations
The aim of this study is to step towards developing a universally applicable workflow for processing long-term geoelectrical monitoring data on slow moving landslides that appreciates an evolving geomorphology. For a reliable geoelectrical model, the practitioner must ensure that electrodes are correctly positioned within the geophysical modelling volume that realises the surface geometry (as is the case for near surface geophysics) to ensure accurate modelling of electrical fields inside the imaging algorithm. The duration of the monitoring data available to this study spans 8.5 years, which to the authors' knowledge, represents one of the longest time series analysis of ERI data within the literature, allowing the issues of data quality control, finite element mesh generation and processing time to be explored. To validate the approach we process, geoelectrical monitoring data for a wellcharacterised site, the Hollin Hill landslide observatory (Chambers et al. 2011;Gunn et al. 2013;Merritt et al. 2013) expanding on a previous study by Uhlemann et al. (2017).
A notable advance in this case is the inclusion of electrode and elevation changes into the finite element mesh used to model resistivities. Robust processing of geoelectrical data is necessary for reliable interpretation of the ERI time series and the conversion of geophysical properties to other parameters such as moisture content (Archie 1947) or soil suction (Crawford and Bryson 2018). We anticipate, elements of the movement modelling methodology presented here could be applicable for future hydrogeophysical investigations of landslides. Henceforth, this paper aims to produce i) an efficient solution for interpolating landslide movements from a sparse grid ii) time-lapse landslide surface and distortion maps of electrode arrays iii) time-lapse 3D ERI volumes which capture distortions to the surface of the slope, and geophysical parameters (electrical resistivity).

Field site: Hollin Hill
Geological setting The Hollin Hill landslide observatory is situated on a south facing~12°slope composed of Lower Jurassic, Lias group, sedimentary rocks (Fig. 1). The succession is dominated by marine mudstones, and the stratigraphy of the field site in ascending order is the Redcar Mudstone (RMF), Staithes Sandstone (SSF), Whitby Mudstone (WMF) and Dogger (DF) Formations. The field site is located in the southern part of the Howardian Hills, North Yorkshire, UK, near to the town of Malton (Fig. 1). The background orthomosaic in Fig. 1 was reconstructed from a fixed-wing Unmanned Aerial Vehicle (UAV) survey in May 2016 as described in Peppa et al. (2019). Here the WMF is the actively failing unit, and is observed to be landslide-prone elsewhere in the UK as the Lias group is geographically widespread (Hobbs et al. 2005). The unit is composed of interbedded siltstones and mudstones, which often host sideritric ironstone nodules towards its base; towards the top horizon of the WMF represents an erosional unconformity (Hobbs et al. 2005). Merritt et al. (2013) provide further details on the geological setting and geomorphological attributes of the site. According to Hungr et al. (2014), the landslide can be classified as a composite, slow clay rotational slide and earth-flow. Many movements have been attributed to translational movements at the SSF-WMF boundary resulting in lobes of reworked mudstone material accumulating downslope of the WMF outcrop ( Fig. 1).

Instrumentation and previous studies
The investigation of slope movements at Hollin Hill began in 2005; in the following years, several geotechnical and geophysical campaigns have taken place in order to better characterise the hillslope (Chambers et al. 2011;Merritt et al. 2013). These efforts culminated in setting up a permanent observatory for studying landslide processes with state-of-the-art instrumentation. The first of these instruments was the Automated Timelapse Electrical Resistivity Tomography (ALERT) instrument (Kuras et al. 2009;Ogilvy et al. 2009) which recoded data from March 2008 (when it was installed), and ran almost continuously up until December 2018. The electrodes were arranged in five parallel lines, 9.5m apart (Fig. 2), with an initial interelectrode spacing of 4.75 m. Each line has 32 electrodes (160 in total) buried at 0.1-m depth to protect the array from animals and general field activities. The 3D monitoring array ( Fig. 1) was set up to characterise the hillslope from head to toe and capture resistivity changes in two flow lobes (Fig. 2), referenced as the eastern and western flow lobe in this study. A grid of wooden marker pegs (45 in total) were installed at the ground surface at 20 m intervals located on the surface above electrode lines (Fig.  1). Alongside the ALERT system, several piezometers, tilt meters, and shape acceleration arrays (SAAs) were installed (Uhlemann et al. 2016a); a weather station which is part of UK COSMOS network (Stanley et al. 2019;Zreda et al. 2012) was installed on the stable part of the slope in 2014 (Fig. 1).

Reactivations
Previous studies and surveys of marker pegs show there have been two major reactivations at Hollin Hill. In November 2012, tilt meters recorded displacements on the western flow lobe which corresponded to an unusually wet summer (Uhlemann et al. 2017), with activity ceasing in February 2013. Additionally, a rotational failure was observed just to the east of the monitored area and captured electrodes on the easternmost part of the array (line 5 in Fig. 2a). Uhlemann et al. (2017) found moisture contents derived from electrical resistivity to be comparatively higher than that recorded for previous years, suggesting that the increased moisture content was driving movements. Over the monitored period the easternmost side of the monitoring array has periodically been reactivated, with lateral displacements up to 8.6 m measured by August 2018.
UAV surveys and passive seismic records (unpublished study) show that another rotational back scarp developed within the monitoring array in late April 2016 at the head of the slope, which spanned four of the array lines (1 to 4). From 2016 onwards the rotational back scarp has continued to grow, presently up to 2.5 m deep, spans 36 m in the Easting direction, and~1 2 m in the northing direction (see Figs. 1 and 2 for backscarp location).

Recording geomorphological changes
Approximately every 2-3 years, terrestrial LiDAR (light detection and ranging) scans and UAV photogrammetry surveys have been conducted in order to capture the changing topography of Hollin Hill. Both techniques are suited to site scale investigation and yield DEMs which can be used to estimate surface changes and generate modelling volumes in ERI. Terrestrial (or ground based) LIDAR is a well-established tool for monitoring rock falls and natural slope movements, be it through permanent monitoring solutions (Lingua et al. 2008) or repeated surveys (Delacourt et al. 2007;Guerin et al. 2021;Palenzuela et al. 2016;Rosser et al. 2007). Recent advances in structure from motion (SfM) photogrammetry have yielded centimetric resolutions such that they are comparable to terrestrial LiDAR scans and both are suited to the purposes of ERI. Recently, Peppa et al. (2019) demonstrated repeated UAV surveys as a means to map landslide movements and geomorphological evolution. In addition, marker pegs were surveyed every 8-12 weeks with real-time kinematic (RTK)-Global Navigation Satellite Systems (GNSS). Over a 10-year period, lateral displacements of up to 8.6 m were recorded on the landslide, whilst vertical displacements up to 2.5 m were observed.

Methodology
Time-lapse ERI processing of the ALERT data is complicated by the dynamic surface topography present at Hollin Hill: electrodes  have moved with the landslide, and their position cannot be directly measured given that they are buried. In addition, metrescale geomorphological features have developed at the site during the monitoring period, the rotational back scarp feature spanning array lines 1 to 4 was not present for previous studies of the site (Uhlemann et al. 2017). From a geoelectrical processing prospective, a robust approach to modelling the geoelectrical measurements should be adopted as discrete changes in topography could mask hydrological changes (as electrical current flow will be incorrectly modelled); in our approach, the 3D surface in the ERI modelling volume (and electrode coordinates) is updated according to the movements of GPS peg markers, the overall workflow is illustrated in Fig. 3 and can be summarised in 3 parts: 1. Update known peg positions after an RTK-GNSS survey, superimpose slope movements on to a reference DEM to create a time-lapse surface. 2. Parse geoelectrical data (from the ALERT instrument) and perform quality analysis, this includes applying appropriate filters to raw data. 3. Combing the outputs of steps 1 and 2 to link the geomorphology of the landslide to the ERI models. Firstly, the time-lapse DEM to inform ERI mesh/modelling volume generation and secondly, the geoelectrical data is inverted to produced volumetric image of the resistivity distribution for a given time step.

Digital terrain models
The DEMs used in this study (Table 1) had already been processed, with the effects of vegetation and other artefacts removed, as part of previous research obtained with a fixed-wing UAV (Peppa et al. 2019) and unpublished studies. All UAV-derived DEMs referenced here are described in Peppa et al. (2019). In brief, aerial images were acquired by a fixed wing Quest 300 UAV (www.ukspacefacilities.stfc.ac.uk), equipped with either a Panasonic Lumix DMC-LX5 or a Sony a6000 compact digital camera. The resulting point clouds were constructed through SfM photogrammetry (processed using PhotoScan, www.agisoft.com) which were used to generate final DEMs with a maximum ground sampling distance of 3 cm as explained in Peppa et al. (2019). The individual UAV-derived point clouds per survey were translated and orientated to a fixed coordinate system (Ordnance Survey Great Britain 36; OSBG36) with the inclusion of surveyed ground control points. Errors due to erroneous co-registration of subsequent UAV surveys were cross-validated with benchmark GNSS observations. Areas of dense vegetation were filtered out from the UAV-derived DEMs and the vertical error of the point cloud used here is estimated to be on average below 5 cm, which is sufficient for the purposes of ERI.
Of the three terrestrial LiDAR scans in this study, two surveys (2008 and 2009) were acquired using a Riegl LPM i800AH, situated on a tripod positioned at the base and halfway up the slope; the raw point cloud was post-processed in RiProfile (www.riegl.com) to remove artefacts associated with vegetation. The most recent LiDAR scan (2018) was acquired with a Leica Pegasus: Backpack Mobile Mapping Solution (Lieca-Geosystems 2019) and involved a continuous walkover survey of the field site, and subsequently processed in Pegasus Manager.
Although satellite-based methods such as Interferometric Synthetic Aperture Radar (InSAR) can have millimetric resolution and have been used successfully in mapping landslide movements (e.g., Booth et al. 2020), in our case satellite techniques were found to be inappropriate due to the lack of permanent scatterers (Ferretti et al. 2001) and poor data availability during the relevant time periods.

Movement modelling
The marker pegs were surveyed using a Leica System 1200 RTK-GPS at a higher temporal frequency than the acquisition of DEMs (shown with vertical lines in Fig. 4), hence an interpolation scheme allows the surface of the DEM to distort with the change in the peg position without the need for frequent DEM surveys (Eq. (1)). Figure 4 (a) shows the frequency of peg and number of transfer resistance measurements passed to the inversion scheme, any missing pegs are assumed to occupy their last known position.
We adopt a thin plate spline approach to map lateral and vertical movements on the hillslope from a discrete set of points (surveying pegs in this case) for each time the points have their position recorded. The displacement for any point within a grid square of the surveying pegs is given as Where d is the displacement vector in the vertical and lateral directions at coordinate (x, y), the other parameters denoted a, b and c are model vectors. This can be solved as a system of linear equations such that Where for example dx 1 is the displacement in the x-direction at position x 1 . Note, Eq. (2) cannot be computed directly as the problem is underdetermined, so Lagrange multipliers are needed to solve the system of linear equations. Thin plate splines are well suited for modelling movement on Hollin Hill as they are valid for an irregular grid (Wahba 1990), making use of four points of reference. Previously Uhlemann et al. (2017) used a piecewise planar approach, which used the three nearest reference points to interpolate the electrode movements. However, this method does not produce smoothly varying displacements across the entire grid squares (Fig. 5). Here, Eq. (2) is solved in order to firstly estimate any electrode positions for a given peg survey, and secondly estimate displacements in the DEM. Some parts of the slope are not subject to movements, as observed from repeated field observations and peg surveys, this is the case for the outcrop of the Staithes Sandstone Formation (Fig. 2). Therefore, an additional constraint is placed on the interpolations of electrode positions, such that electrodes downslope of the flow lobes are fixed (Fig. 2a).
For each time the pegs were surveyed, slope movements are modelled to produce a time series of electrode coordinates and DEMs. When a LIDAR or UAV survey took place during the monitoring period, the reference DEM is updated (Table 1). During the monitoring period, any broken or missing surveying pegs were replaced in-situ, and at no time were the pegs returned to their starting positions (hence the interpolation scheme works on an irregular grid). In order to maintain consistency between DEMs for time-lapse analysis, the point cloud from each UAV or LIDAR scan (Table 1) is filtered with a 2D 1 × 1 m moving average window with 0.5 m tolerance, to avoid interference from vegetation features left inside the DEMs. The point clouds are then downsampled on to a regular 1 m grid for the purposes of ERI mesh generation using a bilinear interpolation scheme. As 3 different field and processing techniques were used to acquire each DEM, the point clouds are then aligned using CloudCompare (GPL-software 2020) against the DEM acquired in July 2008.

Measurements
The ALERT system was set up to record multichannel dipoledipole measurements (Binley and Slater 2020), for both in-line and cross-line (equatorial) configurations. Raw measurements are in the form of transfer resistances (TR): the ratio of a difference in voltage between two electrodes and the current injected in the other two electrodes of a specific four electrode configuration. The  dipole lengths, a, on inline measurements range between 1 and 4 electrode spacings (4.75 to 19 m) with inter-dipole separations, na, where n = 1-8 (Uhlemann et al. 2017). Dipole-dipole equatorial measurements are made on adjacent lines, where a = 9.5 m and n = 0.5, 1.0, 1.5, …, 9.5, 10.0 (since the spacing between adjacent lines is twice the along-line spacing). The ALERT data are stored by date and compiled into a time series of 3D ERI data, with 929 entries in total (Fig. 4) between the 21 st of December 2009 and the 16 th of September 2018. Note that although ALERT was installed in March 2008, automated recording did not start until January 2009 and cross-line measurements were not added to the ALERT scheduling files until December 2009. Generally, measurements were made every 2-3 days, however there are data gaps due to power failures and equipment malfunctions associated with ALERT and its supporting infrastructure. ERI data quality varied widely during the monitoring period, this is largely driven by seasonal changes in contact resistances; which were higher during summer months due to the decreased moisture content of the ground surface (Fig. 4b), resulting in poorer galvanic contact between the electrodes and their surrounding material. Consequently, more data are filtered out during the summers (Fig. 4a). Breakages in the electrode cables (as a result of movements) rendered some electrodes inoperable, also contributing to diminished data quality. The ERI cable on line 5 had to be repaired during the monitoring period due to breakages (Fig. 2a, Fig. 4a), which rendered four electrodes on this line inoperable from February 2015 until October 2016. The measurements are filtered out based on five criteria (two of which refer to a measurement of reciprocity (Tso et al. 2017)): 1. A contact resistance over 5000 Ω, as these measurements are likely to have a high signal to noise ratio given that the environment is relatively conductive (apparent resistivities below 200 Ωm are observed). 2. An approximated apparent resistivity outside the range of 0 and 200 Ωm. Apparent resistivity is computed by multiplying the TR measurement by a geometric factor (Binley and Slater 2020), which differs depending on the array configuration, it's value is valid for homogenous flat ground. In this case, positive geometric factors are anticipated due to the geometry of active electrodes, and measurements with over 200 Ω m corresponded with erroneous TR measurements which resulted in artefacts in the ERI inversions. 3. An amplitude ratio (measure of waveform symmetry) outside the range of 0.85 and 1.15. Beyond this desired range shows that the alternating current signal is asymmetrical. 4. A reciprocal error over 10%, which is often taken as a standard cut off for reliable transfer resistance measurements (e.g., Carrigan et al. 2013). 5. Measurements without reciprocals, as their reliability cannot be assessed.>

Fig. 4 (A) Number of transfer resistance measurements retained after parsing the ALERT data with cross line measurements, including measurements in forward and reverse configurations. Dates of cable repairs and peg surveys are also indicated. Summer and autumn months are greyed. (B) Median contact resistance computed for each day ALERT ERI survey was acquired as ALERT reports a contact resistance for each resistivity measurement
Generally, over 3000 individual measurements are retained for inversion (approximately 68% of viable measurements made in the forward direction).

Inversion workflow
As noted by Uhlemann et al. (2017) (and references therein), incorrect electrode positions within the ERI inversion lead to artefacts in the resulting model. As the ERI surveys have a higher temporal frequency than the peg surveys, a linear interpolation (using time and displacement as input) is used to sample the estimated displacements onto the days which ERI surveys took place (every~3 days). This is deemed appropriate as the landslide is slow-moving and significant movements are captured by the GPS surveys of the marker pegs.
For each time step in the time-lapse ERI, a new mesh with unique topography and electrode nodes is generated; this is necessary to realise vertical and lateral landslide movements in the ERI inversions. The options for time-lapse inversion are therefore limited in this case. Difference inversions (LaBrecque and Yang 2001) do not allow for changing meshes or electrode positions, and time-lapse inversions with moving electrodes have only recently been demonstrated for 2D problems (Loke et al. 2014;Loke et al. 2018). Here we adopt a similar custom workflow to that of previous studies (Uhlemann et al. 2017;Whiteley et al. 2020) through a baseline-constrained approach. A nearest neighbour lookup is used to translate the baseline model values onto each time-lapse mesh. Compared to Uhlemann et al. (2017), who considered a shorter time series, the inclusion of topography is necessary due to the surface changes during the timescales of this study, and whereas a fine mesh was used such that electrodes could move on the same mesh for each ERI inversion, here a coarser mesh is used during the inversion to give a comparatively modest computation time. Of the 929 datasets collected in total, 914 were processed (the rejected surveys have fewer than 500 valid measurements). The inversions were run on a high-performance cluster, across two Intel nodes each with 16 logical processing cores, taking approximately 3 days to run. The baseline inversion is taken from 15 th of April 2010 as this represents a time of intermediate saturation on the hillslope, good data quality and when the landslide was not influenced by movements (Uhlemann et al. 2017).
We use E4D (Johnson et al. 2010) for the 3D ERI on a tetrahedral mesh, as the code scales with computational resources. Additionally, the ResIPy python code is used to prepare data for inversion (Blanchy et al. 2019). Weighting the measured transfer resistances by a reciprocal error model has been shown to produce more robust results (Tso et al. 2017), therefore for each ERI step a unique reciprocal error model is computed based on multi-bin analysis (Binley and Slater 2020;Mwakanyamale et al. 2012). The average resistance of both the forward and reciprocal measurement is taken into the inversion, and the absolute reciprocal error is defined as where R for and R rev are transfer resistance measured in forward and reverse mode, respectively. A different error model is required for each ERI time step because of the different error characteristics present for different seasons and data quality present in the time-lapse data. For all inversions, a constant 2.5% is added to the reciprocal errors to represent the forward modelling errors, as this was found to result in a spatially and temporarily smooth model comparable to previous investigations (e.g., Merritt et al. 2013). Resistivity is expected to vary smoothly from the baseline (Uhlemann et al. 2017), and hence an L2 norm (Loke et al. 2014) is applied as a temporal constraint as well as a smoothness constraint. For time-lapse inversions a relative weight of 0.1 is used as a baseline constraint verses 0.9 for spatial constraint; this encourages a smooth spatial result over smoothed temporal changes. E4D was assigned a target Chi-squared (χ 2 ) value of 1.1. These parameters were found to minimise inversion artefacts, whilst converging on reasonable solutions. We use a custom mesh generation scheme, a flat tetrahedral mesh is generated within Gmsh (Geuzaine and Remacle 2009), and the topography is transposed onto the mesh using triangulation interpolation. Mesh node boundary conditions are then computed, such that the upper surface of the mesh is considered a zero flux boundary (i.e., cannot transmit electrical current) and the mesh is exported into the tetgen format (Si 2015) used by E4D. The baseline inversion is done on a finer inversion mesh to encourage accurate as possible starting resistivities for each subsequent time step during nearest neighbour lookup.

Temperature correction
Electrical resistivity varies as a function of temperature, consequently time-lapse ERI volumes should be corrected for changes in seasonal temperature to avoid misinterpretation of inversion results that could otherwise be confused for hydrological changes (Chambers et al. 2014). The same seasonal depth, z, and temperature model is used to correct the inverted resistivities (post processing) here as in previous studies (Uhlemann et al. 2017) with T mean as the average annual air temperature, T as the difference between the largest and smallest annual temperatures, φ is a phase offset to bring surface and air temperature into phase, d is a characteristic depth and is defined as the depth where T has decreased by 1/e (Brunet et al. 2010). t is the day in the year. The depth of the (barometric) centre of each cell in the mesh is computed, and the corrected resistivity calculated using the ratio model (Ma et al. 2011;Uhlemann et al. 2017) expressed here in terms of resistivity where ρ is the cell resistivity at temperature T model , α is the temperature correction factor, set at −0.02°C −1 , and T ref is a constant reference temperature, in this case 20°C. The constants used in Eqs. (4) and (5)

Landslide kinematics
The modelling of the electrode movements allows for an assessment of landslide kinematics at Hollin Hill (e.g., Hutchinson 1983) over a 10-year period, Fig. 6 illustrates the direction and relative magnitude of lateral electrode movements. At the end of 2012, downslope movements have been observed on the eastern flow lobe, which are accompanied up slope with a rotational failure just to the east of the monitoring array, affecting the electrodes on line 5. These movements correspond to the reactivation of the eastern flow lobe documented by Uhlemann et al. (2017) (Fig. 6b, c and S1). We compute differences from the baseline elevation measured in 2008 (Fig. 6a). The back scarp feature (spanning lines 1 to 4) is clear in the elevation models after April 2016: a decrease in the surface elevation is observed, whilst downslope of the scarp an increase in surface elevation occurs. This supports an interpretation that a rotational slip plane is present at depth; at the head of the failure a slump can be observed corresponding to an accumulation of material. After the development of the rotational back scarp feature (Fig. 6) downslope displacement of the electrodes can be observed, as electrodes move with the rotating mass. Variations in elevation associated with the eastern flow lobe show a decrease near the crest of the lobe, and an increase at the head of the lobe, indicating a downslope translation of material that accumulates at the head of the flow lobe (Fig. 6d). Field observations support this hypothesis, as freshly disturbed material can be observed at the toe of the flow lobe. Although movements in the midsection of the sliding material are slight, it is also likely that this material feeds the flow lobes which move in turn. Any movements downslope on the lobes reduce the support for the upper part of the slope and further encourage the development of rotational back scarps. The observations of electrode and slope movements (as well as field visits) show that failure at the top of the slope is progressing westward in the field area. Although beyond the scope of this study, the electrode displacements effectively provide a 3D displacement field which can be quantitatively accessed to map the slip surface at depth (Aryal et al. 2015;Booth et al. 2020).
Through the workflow described in the foregoing ("Methodology" section and Fig. 3), an ERI time series is produced where features such as the back scarp and flow lobes evolve naturally in the ERI inversion mesh. Figure 7 shows the development of a rotational backscarp at the head of the landslide being reproduced in the time-lapse modelling mesh which captures the inverted resistivities. The depth of the backscarp feature grows from April 2016 to December 2016.

Inverse model validation
The statistical validity of inverse models are generally accessed through a value of χ 2 (Constable et al. 1987;Günther et al. 2006), in the ideal case that there are no modelling errors and data errors are fully realised a χ 2 of 1 should be obtained (Johnson 2014). In this case, E4D converged on a target χ 2 of 1.1 for each time-step showing reasonable fit between ERI models and the ALERT data. Note that setting a target χ 2 of 1 meant E4D could not achieve convergence for all timesteps, it can be expressed as the model misfit over the number of measurements, N, as: where W d is the data weight vector (obtained from the reciprocal error model), d is the measurement vector and f(m) is the forward response to the model parameters m (Binley and Slater 2020). To assess the reliability of results, we ran a separate baseline constrained inversion for March 2017 where: both electrode coordinates and topography are updated, only the electrode coordinates are updated, and neither the topography and electrode positions are updated (Fig. 8), respectively these are referred to as the updated, partially updated and none updated inversions. The percentage RMS (root mean squared) error for these respective inversions are calculated as 6.9, 7.2 and 8.2%, and hence the updated inversion (using the proposed workflow) yielded the best fit in this case. RMS in this case is defined as: where R obs and R sim are the observed and modelled transfer resistances, respectively, and N is the number of measurements. For the partially updated and none updated inversions there is a negative resistivity anomaly present on the eastern flow lobe, which is not consistent with expected resistivity changes or updated inversion. Furthermore, the updated inversion shows positive changes in resistivity compared to the baseline inversion implying relative drying, however the partially updated inversion shows an overall negative change in resistivity implying relative wetting, altering the hydrological interpretation of hill slope processes. This demonstrates the importance of topographic variations when interpreting subtle changes in resistivity as hydrological changes maybe masked if topography is not updated in the time-lapse inversion volumes.
In October 2020, the ALERT system was fully decommissioned and the buried electrode arrays were recovered in preparation for the installation of a new monitoring system. It is challenging to assess the success of the interpolation scheme given the electrodes had been placed in the ground 12 years prior, as landslide movements, particularly in the flow lobes, made it difficult to relocate electrodes which had become disconnected from the buried cable. Additionally, many of the original pegs had perished and hence the quality of the interpolation likely suffered. Where the electrodes were found in place an RTK GPS was used to survey their final position, of the original 160 electrodes 108 (67.5%) were recovered. On average electrodes predicted by the interpolation were 0.51 m away from their final recorded position, the median value is 0.25 m and the standard RMS between the predicted and observed displacements is 0.79 m. For context, the electrodes moved 1.79 m on average and the maximum observed displacement was 7.88 m, suggesting a reasonable fit between the interpolated and observed electrode positions.

Discussion
Time-lapse ERI data is difficult to visualise for multiple time steps in a static format, therefore results of the ERI workflow are presented in S2. To summarise, overall increases in resistivity can be observed during the summer months, which is associated with decreased moisture content. Elevated moisture contents during winter months are associated with lower resistivities. This is in accordance with seasonal moisture content variations observed by Uhlemann et al. (2017). One approach to assess spatial and temporal variability in the resistivity results is to calculate a coefficient of variation (standard deviation of each point in its time series over its mean). Although each mesh in the time series is different, a nearest neighbour lookup scheme can be used to map cell resistivity values onto a representative mesh, from which statistical analysis can be made as in Fig. 9. The Whitby mudstone downslope of the back scarp (in a rotational slump) experiences relatively little change compared to the flow lobes or back scarp area whilst maintaining a relatively low resistivity, indicating the slump retains a high level of moisture throughout an annual cycle. This supports previous interpretations of the hillslope hydrogeology for the Hollin Hill landslide that included perched aquifers (Gunn et al. 2013;Uhlemann et al. 2017). It is likely positive pore pressure under the slump encourages movement on a slip plane at depth. Significant changes in resistivity on the flow lobes (Fig. 9) can be attributed to relative drying during

Workflow
Through monitoring geomorphological changes, it is possible to interpret slope failure mechanisms (Hutchinson 1983); in this case, a rotational failure is observed in slope movements. Consequently, the inclusion of time-lapse DEMs likely improves the quality of inverted images (Fig. 8) as the modelling of the potential field during the imaging process is sensitive to surface topography, for example the development of the backscarp (over lines 1-4) would be particularly troublesome for conventional time-lapse ERI. The approach adopted here facilitates a two phased interpretation through i) visualisation of slope movements characterising the external nature of the landslide though time and ii) capturing the internal structure of the landslide through volumetric electrical imaging which by extension can be related to moisture contents.
With the exception of certain models in the ERI series, which are associated with poor raw data quality (particularly in 2018 when the number of TR measurements drops off significantly), the time-series analysis could be taken further with petrophysical relationships between resistivity, moisture content (Uhlemann et al. 2017) and other critical parameters for assessing slope stability such as soil suction (Crawford and Bryson 2018). Alternatively, more involved workflows could conceivably couple hydrological and geoelectrical modelling through petrophysical relationships (Johnson et al. 2017;Revil et al. 2020), allowing for robust assessments on the slope hydrogeology through time.
A limitation of the workflow proposed here that it fails to account for sudden changes to the slope surface which has been recorded by SAAs and tilt meters, rather treating alterations to the slope surface as smoothly varying between different DEM surveys. We suggest further coupling between in-field sensors, like SAAs, and interpolation of movements to force distortions to slope topography and electrode positions to occur within discrete time windows where movement is recorded. Another drawback is that the approach adopted here relies on repeated field visits that are labour intensive and time consuming, hence an automated approach to monitoring slope movements (Le Breton et al. 2019;Wilkinson et al. 2016) would be beneficial to future studies.

Conclusions
Landslide monitoring through ERI is likely to become more pervasive in coming decades, as the method is suitable for long-term applications and provides volumetric estimations of hydrological parameters that complement more conventional point sensors. However, accurate modelling of the electrical potential field requires a good understanding of the slope geomorphology, which as demonstrated here is subject to slope movements if the landslide is active. Previous papers (Uhlemann et al. 2015;Wilkinson et al. 2016) have addressed modelling electrode movements in ERI inversion, whilst this study proposes a methodology that fully incorporate landslide kinematics (S1) into the inversion workflow.
Time series ERI volumes (S2) capture the changes in slope topography, necessary for avoiding imaging artefacts, and electrical resistivity. Variations in the latter can be reasonably explained by seasonal fluctuations in moisture content observed at Hollin Hill (Uhlemann et al. 2017). Although the 4D ERI data were processed with baseline constrained inversion scheme, reasonable time-lapse results were achieved for the majority of time-steps. Higher χ 2 values associated with inversions where the electrode coordinates or topography are not updated demonstrate that the inversion workflow described here (Fig. 3) improves the quality of inverted results and is necessary for reliable hydrological interpretation of time-lapse ERI volumes; hence, establishing a framework (and corresponding algorithms) for processing (hydro) geophysical datasets resulting from long-term monitoring solutions on active landslides.
Relationships between electrical resistivity and soil moisture are well-documented; hence, geoelectrical model time series can be interpreted in terms of hydro-mechanical parameters through petrophysical calibrations. Furthermore, linking between relevant weather data, petrophysical relationships, ERI data and landslide kinematics (through the framework proposed here) could be used as forcing datasets inside of a coupled hydrological modelling and ERI approach. This is a crucial step forward for developing geoelectrical landslide monitoring techniques and anticipating potential failure events.