Electrical resistivity monitoring of river–groundwater interactions in a Chalk river and neighbouring riparian zone

In the past several decades, there has been considerable interest in groundwater– surface water interactions and their ability to regulate and cycle nutrients and pollu-tants.These interactions are spatially and temporally complex,but electrical resistivity imaging can be a useful tool for their characterization. Here, an electrical resistivity imaging monitoring array was installed laterally across a groundwater-dominated Chalk river and into the adjacent riparian wetland; data were collected over a period of 1 year. Independent inversions of data from the entire transect were performed to account for the changing river stage and river water conductivity. Additionally, data from just the riparian zone were inverted using a temporally constrained inversion, and the correlation between the riparian zone resistivity patterns and river stage was assessed using time-series analysis. The river stage and the Chalk groundwater level followed similar patterns throughout the year, and both exhibited a sharp drop following cutting of in-stream vegetation.For the independent inversions,fixing the river resistivity led to artifacts, which prevented reliable interpretation of dynamics in the riverbed. However, the resistivity structure of the riparian zone coincided well with the intrusively derived boundary between the peat and the gravel present at the field site. Time-series analysis of the inverted riparian zone models permitted identification of seven units with distinct hydrological resistivity dynamics (five zones within the peat and two within the gravel). The resistivity patterns in the gravel were predominantly controlled by up-welling of resistive groundwater and the down-welling of more conductive peat waters following the river vegetation cutting event. In comparison, although the vegetation cutting influenced the resistivity dynamics in the peat zones, the resistivity dynamics were also influenced by precipitation events and increasing pore-water conductivity, likely arising from biological processes. It is evident that such approaches combining electrical resistivity imaging and time-series analysis are useful for understanding the spatial extent and timing of hydrological processes to aid in the better characterization of groundwater-surface water interactions.


I N T RO D U C T I O N
In recent decades, there has been significant interest in the biogeochemical cycling that occurs when groundwater (GW) and surface water (SW) mix in the subsurface (e.g. Hill et al., 1996;Prior and Johnes, 2002;Johnes et al., 2020). GW-SW mixing regulates subsurface temperature and creates unique biogeochemical conditions. Additionally, the presence of mineral surfaces and microbes is such that biogeochemical cycling can actively transfer the nutrients and pollutants across the GW-SW interface, and consequently mitigate ecological degradation throughout the catchment (Kuusemets et al., 2001;Shabaga and Hill, 2010).
However, the importance of biogeochemical transformation may be misrepresented in catchment management literature due to the complex nature and variability of processes occurring in different settings. For instance, the efficiency of biogeochemical cycling is controlled by the sources and timing of water mixing, the substrate and organic matter content, and the microbial communities present, all of which can be drastically different between sites, even within the same catchment (Frei and Peiffer, 2016;Bernard-Jannin et al., 2017). Characterizing these parameters is not trivial, and limited spatial resolution of intrusive methods may make it difficult to unravel complex interactions. There is a need for methods to reveal properties and processes related to GW-SW interactions, for example the timing and extent of GW-SW mixing, with high spatial and temporal resolution.
In recent years, there has been an increasing application of geophysics to characterize GW-SW interactions (McLachlan et al., 2017). For example, electrical resistivity measurements obtained from electrical resistivity imaging (ERI) or electromagnetic induction methods have been used to characterize riverbed structure (e.g. Crook et al., 2008) and to reveal zones of GW up-welling (e.g. Binley et al., 2013). Time-lapse ERI is particularly useful for monitoring GW-SW interactions because changes in resistivity can be interpreted in terms of changing saturation or changing pore-water conductivity. For example, time-lapse ERI has been used with saline tracers to monitor GW-SW interactions over several hours in the hyporheic zone (e.g. Ward et al. 2010;Toran et al., 2012). Additionally, time-lapse ERI has also used natural electrical contrasts (e.g. from solute content or temperature) of water to monitor longer term GW-SW interactions (e.g. Johnson et al., 2012;Wallin et al., 2013;Steelman et al., 2017).
In this work, an ERI array was installed to characterize the timing and location of GW-SW interactions occurring beneath the River Lambourn and the adjacent riparian zone at the Boxford Wetland, West Berkshire, UK. The field site was investigated previously by Uhlemann et al. (2016) who used time-lapse ERI to reveal the complex hydrological patterns involving a peat with two-layer behavior and an underlying gravel. Here, ERI data were collected diurnally over 1 year to monitor electrical resistivity in both the riverbed and the riparian zone, and determine their relation to the changing stage and GW levels. The intention here was to assess the usefulness of ERI monitoring for understanding the timing and spatial extent of GW-SW interactions.

Field site
The River Lambourn and its associated wetlands, comprise one of the least impacted Chalk river systems in the UK, and they have been the subject of many publications to understand groundwater-surface water (GW-SW) interactions (e.g. Allen et al., 2010;Johnes et al., 2020). Furthermore, the Boxford Wetland (Fig. 1) is a Site of Special Scientific Interest (Natural England) and a Special Area of Conservation (EU Habitats Directive) owing to the habitat it provides for flora and fauna; particularly Desmoulin's whorl snail (Vertigo moulinsiana) and river-water crowfoot (Ranunculetum fluitans), which exhibits prolific growth during the spring-summer months. As part of the river management, this aquatic vegetation is cut during the year, primarily to maintain fish habitats for angling and has been shown to significantly influence the river stage and GW levels (Old et al., 2014).
The River Lambourn has a base flow index of 0.96 (Marsh and Hannaford, 2008) at the nearest hydrometric station (6 km SE of the field site); the river stage is therefore predominantly GW controlled. The River Lambourn catchment, and wider Berkshire Downs area, is underlain by extensive, broadly horizontal, Upper Cretaceous Chalk up to 250 m thick in some areas. At the field site, the underlying Chalk is characterized by variable degrees of weathering such that there is a well-developed, largely impermeable, putty layer in its upper sections. The presence of this putty layer is thought to govern the hydrology of the site; for example, regions where this layer is poorly developed form zones of preferential GW discharge (Younger et al., 1989).
The extent of the putty chalk was investigated by Chambers et al. (2014a) using three-dimensional (3D) electrical resistivity imaging (ERI), and locations of GW up-welling within the wetland were revealed by House et al. (2015) via a temperature survey. At the site, the Chalk is overlain by Quaternary alluvial gravel and peat with a typical combined thickness of 3-7 m. The gravels exhibit braided structures, and the basal gravels are characterized by high proportions of reworked chalk, whereas the peat is interpreted to have accumulated during a series of wet and dry periods, as evidenced by the predominance of organic matter sourced from terrestrial and aquatic plants (Newell et al., , 2016. From temporal resistivity patterns, Uhlemann et al. (2016) observed that the peat comprised two hydrologically distinct layers. It was revealed further using intrusive sampling that these layers were separated by a thin layer of clay. Additionally, Musgrave and Binley (2011) used time-lapse ERI over a 1 year period in an adjacent wetland on the east side of the Lambourn, finding that seasonal temperature variation dominated the resistivity signal. They interpreted areas of low temporal variability to be sites where stable GW temperatures subdued seasonal temperature variations within the subsurface.

Site selection
Before the installation of the ERI array, a site of GW upwelling was identified from a riverbed temperature survey. A temperature probe was built by housing a thermocouple in a 1-cm diameter steel pipe, 2 mm holes were drilled in the vicinity of the thermocouple and the end of the pipe was flattened and sharpened to enhance riverbed penetration. The probe was driven into the upper 5-10 cm of the riverbed during winter (08 Jan 2017); throughout the survey, the river water temperature was 7°C and the GW was 10.5°C within the Chalk piezometer ( Fig. 1). Riverbed temperatures ranged from 7 to 10.5°C, and locations with temperatures close to the temperature of the Chalk GW were interpreted to be areas of GW up-welling. The ERI array was installed to coincide with a zone of elevated riverbed temperature (Fig. 2), i.e. a zone of perceived GW-SW connectivity. Additionally, although there is other potential zones of GW-upwelling, this site was chosen beacuse of the relatively low river stage, non-complex bedform geomorphology, and proximity to previously installed instrumentation.

Electrical resistivity imaging data acquisition
ERI data were collected using a PRIME ERI monitoring system (Huntley et al., 2019). The electrode array comprised Locations of electrodes (white dots) and data loggers, view is in the downstream direction. One set of temperature loggers was installed on the East bank (EB-1) and three were installed in the West bank (WB-1, 2 and 3). EB-R and WB-R were used to log stage, temperature and electrical conductivity of the river. WB-G was used to monitor the temperature and water level in the gravel. Logger sets WB-1, 2 and 3 and EB-1 were installed < 0.5 m either side of the ERI transect, EB-R and WB-R were installed ∼4 m downstream and ∼1 m upstream of the ERI transect, respectively. 32 point electrodes each ∼0.55 m apart; a larger spacing was used at either riverbank to prevent electrodes being periodically exposed to the air during low river stage (Fig. 3). Two electrodes were positioned in the east bank, 18 on the riverbed, and 12 in the west bank; concrete slabs were used to prevent movement of the cable throughout the survey period. A dipole-dipole measurement sequence was selected as they provide good vertical and horizontal resolution (Chambers et al., 2002). A sequence with dipole lengths (a) of 1-5 and dipole separations (n) of 1a to 6a was used, resulting in a total of 438 measurements along with a full set of reciprocal measurements.
ERI data were collected diurnally from 16 Nov 2017 to 30 Nov 2018, at both 03:30 and 15:30 GMT. Two river loggers were housed in separate PVC pipes and installed in the river to monitor stage, temperature and electrical conductivity; a series of temperature loggers were installed to correct the ERI inversions for seasonal temperature variations (Fig. 3). Additionally, the GW levels in the Chalk and gravel piezometers (see Fig. 1) were monitored with pressure loggers, Chalk water level was monitored from 11 Feb 2018 and gravel level was logged from the start of the ERI monitoring period, with a logger malfunction leaving a period from 08 Feb 2018 to 27 Jun 2018 absent from the record.

Data quality
Before inversion, data quality and general patterns within the ERI data were assessed. Reciprocal errors were calculated for each quadrupole from the difference between direct and reciprocal measurements; a mean resistance was also calculated to express reciprocal errors as a percentage. Apparent resistivity values were calculated to explore the pattern of data through-out the monitoring period by assuming a flat topography and an electrode spacing of 0.55 m.
Data quality was good with most measurements having a reciprocal error < 0.5%. Measurements with reciprocal errors exceeding 10% or apparent resistivities exceeding 400 ohmm were removed from their respective data sets. To ensure that datasets had similar sensitivity patterns, datasets with less than 416 ohmm (95%) mean resistivity values were discarded. This resulted in a total of 80 datasets out of 739 being discarded; this included 25 day period (between 02 Sep 2018 and 26 Sep 2018) where high contact resistances of electrodes in the riparian zone gave data with high reciprocal error.

Time-lapse electrical resistivity imaging inversion
Time-lapse ERI data can be inverted through performing standard, independent inversion of data collected at different times (e.g. Steelman et al., 2017). However, more elaborate methods have also been proposed, such as inversions where data are constrained to a previous model. This can be done by inverting data with some penalty term for divergence from a prior model (e.g. Oldenborger et al., 2007), or difference inversions whereby the difference between two datasets is used to model deviations from a model obtained from one of those datasets (LaBrecque and Yang, 2001). For the difference inversions, systematic noise, for example arising from poorly surveyed electrode positions, poor galvanic contact, and numerical measurement errors can be removed. Such approaches are useful in cases where systematic errors dominate measurement errors. However, as investigated by Lesparre et al. (2017), special care must be taken when characterizing errors in difference inversions.
When large datasets are concerned, sequential inversions may be time consuming; alternatively, datasets may be inverted simultaneously using some temporal regularization to link datasets (e.g. Johnson et al, 2012;Wallin et al., 2013). In this way, data may be treated in parallel, which can be efficient for datasets collected over long periods or with high temporal resolution. Although these methods require a consistent mesh for each inversion, Wallin et al. (2013) demonstrated how the inclusion of a zone of fine mesh elements could be used to limit model regularization across a dynamic water table and improve interpretation of dynamic GW-SW patterns.

River-based electrical resistivity imaging
It is important to note that additional complications arise in aquatic ERI applications, in comparison to terrestrial applications, due to the abrupt transition between a conductive river and a more resistive riverbed. For example, McLachlan (2020) demonstrated using synthetic modelling that fixing the river resistivity in the inversion significantly improves riverbed characterization. However, the increased sensitivity of measurements to the river, in comparison to the riverbed, is such that erroneous fixing of river resistivity can lead to significant artifacts within the riverbed region of the inverted model. This issue becomes more problematic in time-lapse surveys given that erroneous fixing of river regions in models for sequential datasets could produce significant artificial changes through time, especially if river parameters (i.e. stage and conductivity) change.
Due to these complexities, in this work, for inversion of data from the full ERI transect, a new mesh was generated for each dataset to account for the change in river stage and inversions were conducted independently. Although it could be argued that an approach similar to Wallin et al. (2013) with consistent mesh and dynamic water table could be used, this would require the incorporation of a layer with infinite resistivity to represent the air. Alternatively, an approach involving interpolation of prior inversions onto a new mesh for difference inversions could have been used; however, it was anticipated that, if present, resistivity changes related to GW-SW interactions would be revealed by standard independent inversions. Finite element meshes were generated using Gmsh (Geuzaine and Remacle, 2009). Although the meshes for each dataset were inherently different, they were generated using the same characteristic lengths to obtain similar finite element sizes and thus minimize substantial differences in forward modelling errors.

Error modelling and inversion
As noted, in this work reciprocal measurements were obtained to characterize errors. Reciprocal errors are sensitive to both systematic and random components, for example arising from fluctuating contact between electrodes and soil, and the resultant modification of current pathways (Binley et al., 1995). In this work, an error model is used to assign measurement error to ensure appropriate data weighting during the inversion. Reciprocal error (ε rec ) is proportional to the measured resistance (R) and is often expressed by the following relationship: where a and b are fitting parameters. An envelope fit error model was used to encompass the majority of data (e.g. Slater et al., 2000). Log transformed mean resistances were sorted into bins of equal width, and a linear model was fitted between the log of the sum of mean reciprocal error and twice the standard deviation of the reciprocal error, and the mean resistance. This meant that the error model encompassed 97% of the data.

Inversion
Data were inverted using R2 (Binley, 2019), a robust and mature inversion algorithm that uses the L2 norm (of the parameter space) to minimize the misfit. Through reduced local regularization, R2 permits blocking regularization across specified regions, for example the river-riverbed interface. To begin with each dataset was inverted independently, with a river resistivity fixed to the logged value. Several tests were also carried out whereby regularization between the river and riverbed was blocked, but the inversion was able to modify the resisitivity of the river. However, from experience this results in the river being modelled with large ranges in resistivity, especially in the mesh elements near to the riverbanks. Following inversion of data from the entire ERI transect, data obtained from the riparian zone electrodes were inverted separately using a consistent mesh. Furthermore, to ensure smooth changes in modelled resistivity of the riparian zone, data were constrained to the prior inversion using R2 and inversions were conducted sequentially.

Temperature correction
Given the influence of temperature on resistivity, inverted models were corrected for seasonal temperature variations. Using the procedure outlined by Chambers et al. (2014b), temperature data collected from different depths within the peat and gravel were fitted to the following: where T z,t is the average temperature of day t at depth z, T air is the mean annual temperature of the air, A is the yearly amplitude of the air temperature variation, d is the depth by which the amplitude of the temperature variation reduced by 1/e, ϕ is a phase offset, and ω is the angular frequency (2π /365). In this case, the obtained penetration depth of the temperature signal was 1.375 m. Resistivity models were then corrected using the ratio model (see Hayashi, 2004;Ma et al., 2010): where ρ corrected and ρ model are the corrected resistivity and the resistivity obtained from the inverse model, respectively, c is a correction factor, and T target and T model , is the target temperature and modelled temperature. Each inverse model was corrected to the mean annual temperature, 10.09°C, and the c value used here, −2.95°C −1 , was determined experimentally by Uhlemann et al. (2016) for the Boxford field site. The good fit when using both the peat and gravel temperatures indicated that thermal diffusivities of both materials were similar. It is important to note that the corrections applied here are for laterally averaged seasonal changes in temperature and are not for diurnal temperature variations.

Time-series analysis
To aid with the interpretation of the ERI models obtained from the riparian zone data, cross-correlation metrics were used to assess the relationship between the changes in resistivity of each inversion element and the river stage. Summary statistics such as maximum absolute correlation and time lag to maximum absolute correlation offer a robust method for determining areas of the subsurface that exhibit similar behavior. Such analysis was employed by Johnson et al. (2012) and Wallin et al. (2013) to assess the infiltration of water due to the changing river stage of the River Columbia at the Hanford Nuclear Site, Washington, United States. It was anticipated that the increasing river stage would increase the water content of the peat in the riparian zone, that is either from lateral infiltration from the river or vertical upwelling of GW, given the hydraulic connection between the river and GW. Consequently, a negative correlation between river stage and electrical resistivity in the peat was anticipated. However, it is also important to note that if conductive waters are replaced by more resistive waters, a positive correlation could be expected. In addition to obtaining maximum absolute correlations, which indicate the connection of the subsurface to the changing river stage, the lag time to maximum absolute correlation should be related to pore water velocities (Johnson et al., 2012). For example, areas with a short lag time could be used to indicate areas of higher pore water velocity and preferential flow.
Given that there is a gap in ERI coverage between 02 Sep 2018 and 26 Sep 2018, and that correlation analysis requires equally spaced data, only data from 16 Nov 2017 to 02 Sep 2018 was considered. Other gaps in the monitoring (≤ 2 days) were accounted for by linearly interpolating resistivity values of elements from temporally neighbouring inverse models. The cross-correlation between river stage and bulk resistivity of each element was calculated using the Pearson correlation for different time lags, for example: where N t is the number of stage and resistivity values in the time sequence, N k = N t-1 , ρ k is the bulk resistivity at time k,ρ is the mean resistivity of the time sequence, σ ρ is the standard deviation of the resistivity time sequence, S k ,S, and σ S are the corresponding river stage time-series, mean, and standard deviation.

General electrical resistivity imaging data patterns
Over the monitoring period, the river stage ranged from 90.7 to 91.1 m above sea level and rose steadily from mid Dec 2018 due to increasing water level in the Chalk and growth of river water crowfoot (see Fig. 4a). The abrupt drop in stage on 20 Jun 2018 coincides with the removal of the river vegetation; furthermore, this drop was also observed in the Chalk water level (Fig. 4b). The specific conductivity of the river water (∼55 mS/m) was stable throughout the year and matched the specific conductivity of the Chalk groundwater (GW). As a result, the electrical conductivity variation is predominantly dependent on river temperature variation (Fig. 4c). The mean apparent resistivity coincides with the river stage and drops abruptly following the aquatic vegetation cutting (Fig. 4d). Curiously, the aquatic vegetation cutting event results in different patterns of the maximum and minimum resistivities (Fig. 4c): the maximum resistivity shows an abrupt increase following the vegetation cutting, whereas the minimum resistivity shows an abrupt decrease. The increase in maximum resistivity could be attributed to an increase in the measured resistivity of the largest quadrupole spacing given that the sensitivity patterns will be shifted significantly; the fall in minimum resistivity could be attributed to the bulk increase in the conductivity of the river, for example without the presence of plants, and associated trapped sediment (Old et al., 2014). Moreover, it can be seen that the mean reciprocal error also increased during the time immediately following the vegetation cutting. This was due to the increase in contact resistances of the riparian zone electrodes immediately following vegetation cutting, perhaps due to a reduction of moisture content in the uppermost riparian zone.

Figure 4
Seasonal patterns of river properties and ERI measurements: (a) river stage, (b) Chalk water table, (c) resistivity of river water, (d) mean (red points), maximum and minimum apparent resistivity of each dataset and reciprocal error as a percentage (blue points). Water levels are expressed as meters above sea level (masl).

Background resistivity
The resistivity model for the background dataset (from 16 Nov 2017) is displayed in Fig. 5. A two-layer structure can be observed in the right side of the resistivity model, which coincides well with the intrusively derived depths to the peatgravel interface (white dashed line). The bulk resistivity of the peat layer is in the order of 15-50 ohmm, which is in agreement with previous electrical resistivity imaging (ERI) studies conducted at the field site (e.g. Chambers et al., 2014a;Uhlemann et al., 2016). Similarly, resistivities of 70-200 ohmm for the gravels beneath the riparian zone and in the majority of the riverbed also agree with Chambers et al. (2014a) and Uhlemann et al. (2016). However, immediately beneath the river, modelled resistivities exceed 1000 ohmm, for example beneath electrode 5 and between electrodes 19 and 20 in Fig.  5 (note that the colour scale is limited to 1000 ohmm). These extreme values are similar to those investigated by McLachlan (2020) and are likely to have arisen from a combination of reduced sensitivity in the upper portion of the riverbed and complications in fixing of river parameters in the inversion.

Temporal resistivity patterns
The extreme resistivity values beneath the riverbed are particularly noticeable when comparing data collected at two different times. For instance, comparing the background dataset (16 Nov 2017) and the dataset collected the day before river vegetation cutting (19 Jun 2018), changes exceeding ± 150% can be observed in the riverbed (Fig. 6b). Given the comparable specific conductivities of the GW and the river water, and the fact that changes due to temperature could only account for ∼3% change in resistivity per°C, these patterns are likely to be artifacts. As subtle changes in the riverbed resistivity dynamics are likely to be obscured by these inversion artifacts, the focus of the remainder of the paper is on characterizing the riparian zone dynamics.

Time-series analysis
Time-series analysis of the river stage and bulk resistivity for each mesh element of the riparian zone data was used to identify areas that exhibited similar behaviors. Plots of maximum absolute correlation with river stage and time to maximum absolute correlation with river stage are presented in Fig. 7, and several time-series from elements of interest are presented in Fig. 8. The mean resistivity and coefficient of variation for the resistivity for each mesh element throughout the monitoring period were also calculated. The mean resistivity (Fig. 7a) of the entire monitoring period shows the same pattern as [Correction made on 15 July 2020, after first online publication: Figure 6 was previously a repetition of Figure 5 and has now been corrected in this version.] i.e. a more conductive peat overlying a more resistive gravel layer. The highest variation in resistivity can be found in the upper portion of the peat; however, the gravel resistivity also demonstrates some variability (Fig. 7b). It is also important to note that although correlation analysis was done using river stage data, it can be seen from Fig. 4(a) and 4(b) that the Chalk water level and river stage are coincident with one another; the river stage was used in the time-series analysis due to its longer period of coverage. It is immediately apparent that the maximum absolute correlations of peat resistivity and river stage exhibit a negative correlation, whereas the maximum absolute correlation of the resistivity in the gravels is positively correlated with the river stage (Fig. 7c). Furthermore, in the gravels, the most common lag to maximum correlation was 0 days (Fig. 7d) indicating the behavior of the gravel resistivity and river stage coincide well. However, a zone where lag time to maximum correlation is in the order of 30-40 days can be seen at the horizontal position of ∼14.8 m. In comparison, the peat lag times are more disparate, with two prominent zones characterized by short lag times (<5 days).
Given the anticipated high permeability of gravel and the presence of the zone characterized by lag times of 30-40 days, it is important to explore what these cross-correlation metrics actually show in terms of the resistivity and stage timeseries (see Fig. 8). For instance, element 1 (Fig. 8a) exhibits a negative correlation (−0.812) and a large lag time to minimum correlation (66.5 days); it is evident in comparing the time series that although the stage drop following vegetation cutting coincides with a resistivity drop, the river stage does not coincide with the pattern of the resistivity well. In comparison, element 2 (Fig. 8b) is characterized by a negative correlation and a shorter lag time to the minimum correlation (6.5 days). However, it is evident again that the resistivity time-series is not simply shifted as the maximum stage and minimum resistivity occur at similar times. Similarly, complementary bulk resistivity patterns are also evident for elements 3 and 4, but the time lag to maximum absolute corre-lation is substantially longer than the lag between the maximum stage and the minimum resistivity ( Fig. 8c and 8d). Furthermore, although the stage drop following the river vegetation cutting coincides with an increasing resistivity in element 4, the change resistivity pattern during Nov 2017 to Mar 2018 behaves differently from the stage, indicating some other factors may be governing resistivity patterns in the peat (Fig. 8d).
For the elements in the gravels, element 6 is characterized by a positive correlation (0.976) and a lag of 0 days and follows the same pattern as the river stage (Fig. 8f). In comparison, element 7 is characterized by a similar positive correlation (0.91), but a long lag to maximum correlation duration (Fig. 8g). On inspection of Fig. 8(g), it can be seen that following an initial drop in resistivity immediately after the vegetation cutting the resistivity increases above the precutting value. This, again, gives evidence of additional controls on resistivity occurring in the riparian zone, other than the changing river stage. Lastly, although element 8 is located in the gravel, according to intrusive measurements, it behaves more similarly to mesh elements in the peat. This could be because intrusive measurements were not made directly along the ERI line, due to the presence of tree roots, or due to some additional three-dimensional (3D) effects influencing the inversion.
For the majority of the peat (i.e. areas showing a negative correlation in Fig. 7c), the minimum resistivity occurs immediately before the vegetation cutting (Fig. 8). Furthermore, whereas the stage decrease is abrupt, the increasing resistivity following the vegetation cut is more gradual. This implies that the behavior of the stage and resistivity of the peat are not linearly related. Additionally, the metric of time-lag to maximum absolute correlation is somewhat misleading as it can be seen clearly in Fig. 8(b) and 8(g) that it does not give a true measure of the lag time between the river stage and bulk resistivity patterns. It does still, nonetheless, offer a useful metric of identifying areas in the riparian zone with similar behaviors.

Resistivity patterns in the riparian zone
Based on the cross-correlation metrics discussed above, the subsurface was split into seven zones, five within the peat (P1-5) and two within the gravel (G1-2). The thresholds used to discriminate between different zones are displayed in Table 1 and the respective zones are shown in Fig 9(a). The ge-ometric mean resistivities for each zone throughout the monitoring period were then calculated and are shown in Fig 9(b) alongside river stage and precipitation data.
In Fig 9(b), it can be seen that the P1 resistivity gradually increased from Nov 2017 to its maximum resistivity during Feb 2018, following this, resistivity began to decrease. The resistivity of P1 responds to the vegetation cutting by exhibiting a sharp increase, then from Jun 2018 to Nov 2018 it increases Table 1 Summary of the different zones of the riparian zone as defined by cross-correlation metrics and coefficients of variation. P zones denote zones within the peat and G zones denote zones within the gravels.
steadily with some rapid drops followed by gradual increases, which coincide with precipitation events (Fig. 9b), for example on the 14 Oct 2018 and 10 Nov 2018. In comparison, the resistivity of P2 increases slightly from Nov 2017 through to Feb 2017 where it begins to decrease to a minimum around the time of the vegetation cutting. The resistivity of P2 is also perturbed by the vegetation cutting but is characterized by a small drop in resistivity, before returning to pre-cutting values immediately after the vegetation cutting. Following this, the resistivity increases, similarly to P2 and it also experiences sudden drops in resistivity on 14 Oct 2018 and 10 Nov 2018, albeit the signal is more subdued. The resistivity of P3 follows a similar pattern as P2 but it increases just before the vegetation cutting, and immediately decreases following the vegetation cutting, before resuming similar patterns as those in P2. The bulk resistivity patterns for P4 are the most variable; however, they become more stable as resistivity falls, probably due to increasing saturation; in general, the resistivity patterns are similar to those in P1. P5 shows a much less variable resistivity pattern; however, it is still affected by the vegetation cutting and shows a slight, and gradual, increase following the cutting. Also, P1, P2, P3 and P4 all show a sharp increase on 01 Mar 2018; this does not appear to coincide with precipitation; however, it does precede a sudden increase in river stage and could instead be related to changing groundwater (GW) levels.
The resistivity of G1 follows very closely the pattern of the rising stage and is characterized by a drop in resistivity immediately following the vegetation cutting, similarly to P3. From Jun 2018 onward, the resistivity of the zone remains relatively stable, with a few perturbations around early to mid Oct 2018. Similar resistivity patterns are seen in G2; however, the resistivity is both lower and the ranges in resistivity are not as great. Also, immediately following the vegetation cutting the resistivity of G2 increases to a value above the pre-vegetation-cutting level; this is not seen in G1. to increasing pore-water conductivity, perhaps as a result of vegetation and microbe related processes. The increasing porewater conductivity in P1 also agrees with the patterns in P2, G1 and G2 where their bulk resistivity decreases immediately following the vegetation cutting, as a result of an abrupt draining of more conductive water from P1. The resistivity patterns in G1 and G2 following the vegetation cutting are also distinctive. While the resistivity of G1 remains stable following its initial drop the resistivity of G2 increases perhaps due to influx of more resistive water, perhaps as GW levels equilibrate following vegetation cutting. It is evident from Fig 9(b) that the vegetation cutting dominates the resistivity dynamics of the riparian zone; however, patterns within the P1, P2, P3 and P4 indicate that precipi-tation events influence the saturation of the peat layers when resistivity is high, i.e. when moisture content is low. For instance, the influence of precipitation is most evident during Nov 2017 to Mar 2018 and Oct 2018 to Nov 2018. It is evident from Fig. 10 that rainwater infiltration is focused on the top-left portion of the riparian zone following the precipitation event that occurred between electrical resistivity imaging (ERI) measurements made on 13 Oct 2018 at 15:30 GMT and 14 Oct 2018, also at 15:30 GMT.

Determining hydrological properties of different zones
In addition to interpreting the patterns observed in each zone, the information about the hydrological properties of each can also be inferred. For instance, P1 and P4 are seemingly closely connected given that they respond similarly to changes in the river stage and precipitation events. However, given the higher resistivity and more extreme changes in response to precipitation events, it is anticipated that P4 is characterized by higher porosity. Furthermore, given that the resistivity of P4 is always higher both zones are unlikely to have ever been completely saturated during the monitoring period.
In comparison, the resistivity of P5 remains very stable during the monitoring period, and whilst the effect of the vegetation cutting is seen on the resistivity, the response is both subtle and prolonged (Fig 9b). This implies that this zone is not well-connected to the rest of the peat, and pore-water velocities are low. This could be attributed to a higher clay content which would be in agreement with its low resistivity. P2 and P3 both show a decrease in resistivity following the vegetation cutting (i.e. the opposite behavior of P1 and P4). This implies that although throughout the majority of the year there is a limited connection of P1 and P4 with P2 and P3, when hydraulic gradients are large enough, water can be exchanged between zones. The increase in resistivity in P3 immediately before the vegetation cutting is also worthy of mentioning and could be related to the up-welling or lateral infiltration of more restive water after some threshold water level is surpassed.
Lastly, the different behaviors of the gravel zones G1 and G2 are important to discuss given that they were anticipated to exhibit similar patterns. The behavior of G1 is intuitive; the resistivity increases due to the up-welling of resistive GW, it then drops following the vegetation cutting due to the draining of more conductive pore waters before its resistivity stabilizes as the conductive pore waters from the upper peat layers are diluted with resistive GW. In comparison, although G2 follows the same pattern it is less resistive than G1 and the resistivity increase during Mar 2018 to Jun 2018 is not as large. Furthermore, whereas the resistivity drop following the vegetation cutting in G1 is abrupt, it is more gradual in G2. This coupled with the increasing resistivity following the initial resistivity drop following vegetation cutting could be a result of slower influx of more conductive waters from above (e.g. due to the presence of P5) and then the up-welling of more resistive GW however the exact mechanisms are unclear.
The presence of different zones within the peat may be attributed to the complex depositional patterns at the field site. For instance, Newell et al. (2015) observed inter-bedded sands within the peat that they propose could be indicative of former channels or over bank deposits. Similarly, using ground penetrating radar data Newell et al. (2015) demonstrated that the upper portion of the gravels is highly structured with clearly defined channel form features that could give rise to the distinct behavior present within the gravels.

C O N C L U S I O N S
Time-lapse electrical resistivity imaging (ERI) was carried out to reveal vertical and lateral exchange of a river with the groundwater (GW). Despite measures to appropriately account for the river stage and resistivity, resolving processes in the riverbed was not possible. Furthermore, the localized nature of some artifacts could arise from poorly surveyed electrodes. It was demonstrated by McLachlan (2020) that accounting for river stage in inversions can be contentious even in relatively straightforward aquatic cases; additionally here the presence of vegetation is also likely to influence the bulk resistivity of the river. For instance, the presence of vegetation and the associated trapped sediment is such that the bulk resistivity of the river is likely lower than the logged river water resistivity, such an effect could also explain the drop in minimum resistivity immediately following river vegetation cutting. Future work ought to be aware of issues surrounding regularization across the riverbed and the bulk resistivity of rivers, including the presence of river vegetation.
Cross-corelation metrics proved useful for simplifying the large dataset and identifying patterns and regions of the subsurface that behaved similarly. The correlation between the resistivity each riparian zone element, and the river stage was assessed; however, it is important to note that the river stage and Chalk GW had coincident patterns. Although it was apparent that there were additional controls of subsurface resistivity, especially in the peat, and non-linear relationships between bulk resistivity and river stage were present, cross-correlation metrics permitted grouping of distinctive hydrological units. Such identification of distinct hydrological behavior in the subsurface could provide a useful tool for locating areas of enhanced residence times or increased mixing, which may have important implications for catchment health. Moreover, although time lag to maximum absolute correlation was shown to differ from the time lag between maximum stage and maximum/minimum resistivity in some cases, it still provided a useful metric to group zones of similar behavior.
It has been demonstrated that ERI monitoring, coupled with time-series analysis, provides a useful method for assessing the extent and timing of groundwater-surface water (GW-SW) interactions within the riparian zone. Future work could rely on a more extensive installation of GW level and bulk conductivity loggers, which could help to interpret patterns of electrical resistivity at a smaller scale and validate some of the observations. Also, approaches aiming to characterize microbiological processes, alongside similar ERI installations, could help to quantify implications of GW-SW interactions for catchment-scale health.

AC K N OW L E D G E M E N T S
This work is supported by the NERC Envision Doctoral Training Program GA/15S/004 S301). The contributions of Jonathan Chambers and James Sorensen are published with the permission of the Executive Director, British Geological Survey (UKRI-NERC).

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author, PM.