Removing the effects of temperature on electrical resistivity tomography data collected in partially frozen ground: Limitations and considerations for field applications

Time‐lapse electrical resistivity tomography (ERT) can be used to image a wide variety of subsurface processes. It is often necessary to remove the effects of seasonal near‐surface temperature variability in order to interpret the signal of interest. Here, we use a synthetic modeling approach to explore the challenges related to removing temperature effects from ERT data in seasonally frozen ground. Electrical resistivity tomography data collection and processing methods that are often used to improve resolution in frozen ground do not appreciably improve the accuracy of temperature corrected data. A sensitivity analysis showed that errors in the temperature corrected data are primarily caused by error in the inverted resistivity models, and that this sensitivity to model error is much higher in partially frozen ground than in unfrozen ground.


INTRODUCTION
Electrical resistivity tomography (ERT) is a geophysical method that uses electrical resistance measurements to estimate subsurface resistivity. Resistivity is sensitive to several factors, including rock type, porosity, water saturation, pore water chemistry, and temperature, making ERT a valuable tool for a wide range of environmental problems (Samouëlian et al., 2005). Additionally, multiple ERT datasets can be collected over time to image dynamic processes-for example, to monitor snowmelt infiltration (French & Binley, 2004), permafrost evolution (Hilbich et al., 2008), coastal aquifer water quality (Ogilvy et al., 2009), contaminant remediation (Wilkinson et al., 2010), CO 2 sequestration (Carrigan et al., 2013), and many others.
One challenge with time-lapse ERT surveys is compensating for the effects of subsurface temperature. Because resistivity is temperature dependent, subsurface temperature variations can interfere with the signal of interest, particularly in time-lapse ERT studies where the goal is to detect relatively small changes in resistivity. In many cases, temperature effects need to be removed, and failure to do so can lead to a fundamental misinterpretation of subsurface processes (Hayley et al., 2009). Several studies have corrected resistivity data for temperature effects in the 0-35˚C range (Chambers et al., 2014;Farzamian et al., 2015;Jayawickreme et al., 2010). Other studies have incorporated temperature corrections under subzero temperature conditions to estimate unfrozen water content (Dafflon et al., 2016;Hauck, 2002). However, there is a knowledge gap when it comes to temperature corrections in partially frozen ground. Compensating for temperature effects in partially frozen ground is complicated by the large increase in resistivity when the ground freezes (Herring et al., 2019). In this paper, we create resistivity models that include temperature effects, simulate ERT data, and assess the feasibility of removing the influence of temperature in seasonally frozen ground. Our goals are (a) to quantify the error in temperature corrected ERT data, (b) to assess sensitivity to different ERT data acquisition and processing strategies that are commonly used to improve resolution of sharp contrasts, and (c) to identify the root cause of errors in temperature corrected ERT data.

Bulk resistivity-temperature relationship
Archie's equation (Archie, 1941) is an empirical relationship used to describe how bulk resistivity ρ bulk (Ω m) changes as a function of porosity ϕ (m 3 m −3 ), water saturation w (m 3 m −3 ), and fluid resistivity ρ f (Ω m), using unitless fitting parameters (the cementation exponent) and (the saturation exponent). A slightly modified version of Archie's equation (Winsauer et al., 1952) is most commonly used in practice and includes another fitting parameter , often called the tortuosity factor: Although the inclusion of fitting parameter is theoretically unjustified, as it creates an inequality when ϕ approaches 1, its inclusion often produces a better fit to observed data by compensating for systemic errors in porosity, fluid resistivity, and temperature measurements (Glover, 2016). Notably, Archie's equation assumes that electrical conduction occurs in the pore fluid and is therefore only valid for clay-free materials. Modified versions of Archie's equation have explored the influence of electrical conduction along mineral surfaces (Revil & Glover, 1998;Waxman & Smits, 2003) but will not be considered here.
Temperature effects can be incorporated into Archie's equation by modifying a few key terms. The remainder of this section will describe the strategy used by Herring et al. (2019) to introduce temperature dependent terms into Archie's equation above and below 0˚C. Above 0˚C, temperature affects the fluid resistivity by changing its viscosity and ion mobility. Fluid resistivity decreases predictably with increasing temperature (˚C): where ρ f_ref (Ω m) is the fluid resistivity at reference temperature ref (˚C), typically 25˚C, and (˚C −1 ) is a constant tem-

Core Ideas
• The temperature signal in time-lapse data can interfere with the signal of interest. • When the ground surface freezes in winter, the resistivity perturbation is large. • Standard data processing methods cannot remove effects of a frozen surface layer. • Temperature correction is ultimately limited by accuracy of inverted models.
perature compensation factor that is typically close to 0.02 for most ionic solutions (Hayashi, 2004, and references therein). At temperatures below 0˚C, two key terms in Archie's equation are affected by temperature: the water saturation and the fluid resistivity. As the ground freezes, pore water turns into ice, effectively reducing the liquid water saturation. However, not all the pore fluid freezes at a uniform temperature. Some water remains liquid even at very low temperatures (<−50˚C) (Yoshikawa & Overduin, 2005) as a thin film along grains (Cahn et al., 1992). The distribution of water during the freezing process is analogous to soil drying (Black & Tice, 1989;Flerchinger et al., 2006;Spaans & Baker, 1996), and the increase in resistivity at subzero temperatures is often modeled as a reduction in liquid water saturation (Dafflon et al., 2017;Hauck, 2002). A temperaturedependent liquid water saturation term wL ( ) (m 3 m −3 ) is therefore introduced that is related to the initial water saturation w0 (m 3 m −3 ) via a unitless relative saturation term r : During freezing, ion exclusion causes solutes to concentrate in the remaining unfrozen pore fluid, which decreases fluid resistivity. Under the assumption that the amount of solute present in solid ice is negligible (Weeks & Ackley, 1986), the decrease in liquid water content is accompanied by a proportional increase in total dissolved solids (TDS) concentration in the remaining liquid water (g L −1 ) compared with the initial TDS concentration 0 (g L −1 ): The relationship between TDS concentration and fluid resistivity at a standard temperature is expressed as where correlation factor [(Ω m) −1 (g L −1 ) −1 ] typically varies between 0.13 and 0.18 (Lloyd & Heathcote, 1985). Combining Equations 1-5, Archie's equation can be modified to account for subzero temperature dependence by introducing (a) temperature-dependent liquid water saturation wL ( ), and (b) a saturation dependent fluid resistivity term (Herring et al., 2019). The general expression for temperature dependent bulk resistivity becomes

Soil freezing characteristic curve
The bulk resistivity of frozen soils depends on the amount of water that remains unfrozen at a given temperature. The soil freezing characteristic curve (SFCC) describes this relationship between liquid water saturation and subzero temperature. A number of different functions have been used to represent SFCCs including power law (Anderson & Tice, 1972), exponential (Michalowski & Zhu, 2006), piecewise linear (McKenzie et al., 2007), and nonlinear piecewise functions (Kozlowski, 2007). An exponential function typically produces a good fit to observed SFCCs for a variety of soil types (Ren et al., 2017) and is commonly used in practice (Ge et al., 2011;McKenzie et al., 2007). This exponential SFCC is described by where w_res (m 3 m −3 ) is the residual water saturation that persists even at very low temperatures, w (m 3 m −3 ) is the water saturation under unfrozen conditions, and is a unitless fitting parameter.

ERT data
The ERT data is measured as a voltage difference between two potential electrodes or as a resistance (Ω), calculated using Ohm's law = ∕ since the input current value is known. However, since resistance values depend on the measurement geometry, it is usually more practical to display data as apparent resistivities. An apparent resistivity ρ a (Ω m) is calculated by multiplying the observed resistance by a geometric factor to account for survey geometry (Telford et al., 1990). Apparent resistivity data are usually then plotted as a pseudosection, where the lateral position of the data point is the midpoint of the quadripole and the relative pseudodepth is a function of electrode spacing and geometry. Following the method of Cockett et al. (2015), the pseudodepths for twodimensional (2D) dipole-dipole and Wenner surveys can be expressed as (8) where AB midx is the midpoint between the two source electrodes and along the x axis, is the location of source electrode along the x axis, and so on. Although pseudosections can give an overall idea of the main resistivity trends in the subsurface, the data must be inverted to get an approximation of the true earth structure.

ERT inversion
Electrical resistivity tomography inversion is the process of using a series of resistance measurements to estimate subsurface resistivity distribution. The standard ERT inverse problem is ill-posed and therefore requires regularization. The goal of the inversion process is to find a model estimate that reproduces the observed data to within the noise levels of the data and is also geologically reasonable. Tikhonov style regularization (Tikhonov & Arsenin, 1977) accomplishes this by minimizing an objective function Φ( ) that usually uses an L2 norm to calculate data misfit (Φ d ) and model misfit (Φ m ) terms (Farquharson et al., 2002): In the objective function, β is a regularization parameter that adjusts the relative weights of the data term and model term, d is a diagonal data weighting matrix, and m is a combination of model weighting matrices. For a typical twodimensional inversion (i.e., a cross-section), m has the form where s enforces smallness (i.e., deviation from a reference model), and enforce first-order model smoothness in the x and z directions, and λ , λ , and λ , are the respective weighting coefficients. These regularization terms can be used to incorporate information about the geologic structures in the inversion process if they are known a priori. For example, if the subsurface has strong horizontal bedding, we may choose to increase the lateral smoothness coefficient λ to produce an estimated model that is more continuous in the x direction. If the locations of boundaries are known (e.g., a horizontal interface, we can multiply λ by a diagonal matrix of weighting coefficients with zeros at the expected boundary locations to avoid penalizing a sharp gradient across the interface (Elwaseif & Slater, 2012).
Linearizing Equation 9 about the reference model and taking the derivative and setting it to zero yields a system of equations where it is possible to solve for a model update Δ : where is the Jacobian matrix of partial derivatives of the data with respect to the model. The model is updated iteratively until the solution converges to meet a stopping criterion. For more detail on implementation, the reader is referred to Pidlisecky et al. (2007).
Other metrics of misfit can be useful when ERT is used to image frozen regions in the subsurface. Rather than using an L2 norm to calculate data and model misfit terms in the objective function, many studies use an approximation of the L1 norm instead (Hilbich et al., 2011;Kneisel et al., 2014;Scapozza & Laigre, 2014). It is common to employ an L1 model misfit norm in a "blocky" inversion because it allows for better recovery of sharp boundaries (Hilbich et al., 2009;Loke et al., 2003). This is particularly helpful for imaging sharp, high-contrast boundaries between frozen and unfrozen ground. As described in Farquharson and Oldenburg (1998), Equation 11 can be modified to enforce other measures of model misfit by adding diagonal data and model reweighting matrix : If an L1 model norm is desired, the diagonal elements of are where and ε is a small value used in the Ekblom norm (Ekblom, 1973) to ensure differentiability. This reweighting matrix must be calculated iteratively; hence, this method is termed iteratively reweighted least-squares (IRLS). Details on the computation of reweighting matrices can be found in Farquharson and Oldenburg (1998) and Sun and Li (2014).

Temperature correction
Removing the effects of temperature on ERT images involves scaling resistivity values to some reference temperature. For data collected in the 0-25˚C range, Hayley et al. (2007) outlined a method to compensate ERT images to a reference temperature based on the relationship in Equation 2: where ρ bulk_ (Ω m) is in situ bulk resistivity at temperature (˚C) that is above 0˚C and ρ bulk_std (Ω m) is the bulk resistivity at the reference temperature std (˚C). std is usually defined as 25˚C (Brunet et al., 2010), as some representative mean ground temperature (Coscia et al., 2012), or as the temperature at a depth beyond the influence of seasonal temperature changes (Hayley et al., 2009).
A general temperature correction that accounts for freezing can be found by rearranging the rock physics relationship described in Equation 6 to solve for the temperature corrected bulk resistivity ρ bulk_ref (Ω m) at a reference temperature of In many cases, it is not advisable to apply a temperature correction directly to the estimated resistivity model. If large resistivity contrasts exist in the subsurface, regularization blurs those features in the estimated model, causing high resistivities to be underestimated and low resistivities to be overestimated. Since magnitudes are not recovered, applying the above petrophysical relationship will result in overcorrection. Hayley et al. (2010) mitigated this issue by applying a temperature correction to ERT data, rather than directly to the model, and found this method reduced artifacts associated with regularization when the subsurface temperature was in the ∼3-18˚C range. Their data temperature correction uses a first-order Taylor series approximation to estimate the temperature corrected data TC obs , under the assumption that the resistivity perturbation (Δ ) that corrects the estimated resistivity F I G U R E 1 Cross-section of the conceptual model of a leaking frack pond. A leaked contaminant plume travels horizontally from the pond towards the plane of the electrical resistivity tomography (ERT) survey model to a reference temperature is small: where obs and est are the observed and estimated data and TC est is the data associated with the temperature corrected model. In this method, the initial inversion of the observed data yields model est and associated data est , applying Equation 16 to est yields temperature-corrected model TC est , and forward modeling TC est produces data TC est .

TESTING SCENARIOS
The data temperature correction from Hayley et al. (2010) has successfully been applied under unfrozen conditions but has not yet been used in studies of partially frozen ground. In this study, we aim to establish the feasibility of applying a temperature correction in regions where the ground freezes seasonally. Although there are many cases where it is desirable to remove the influence of temperature on an ERT dataset, this paper will focus on imaging a conductive plume in seasonally frozen ground. In the hypothetical situation considered here, a conductive plume originated from a lined surface pond that stored highly saline flowback water from fracking. A permanently installed ERT survey line, 40 m from the pond and oriented perpendicular to the direction of groundwater flow, collected data at monthly intervals to monitor the frack pond for leakage. The hypothetical leak occurred in June and horizontally migrated 40 m during the subsequent 1-yr monitoring period, ending with its center of mass directly below the ERT line at the midpoint of the survey. The contaminant plume in this case is spherical, with a diameter of 14 m, and is assumed to move strictly by advective transport (i.e., no diffusion or dispersion) horizontally away from the pond towards the ERT survey line. A schematic diagram of this conceptual model is shown in Figure 1. Although this configuration is obviously simplified, it represents a valid way to test the temperature cor-rection methodology without introducing unnecessary complexity.
To model the contaminant plume migration, we created time-varying three-dimensional (3D) resistivity models that included the migrating contaminant plume and the effects of temperature. We used numerical modeling to simulate the data measured by a permanently installed ERT array. Before testing whether the temperature correction worked on data collected in seasonally frozen ground, we chose to validate the temperature correction on another scenario where the ground never froze (Scenario A). Therefore, we present two scenarios here: Scenario A, in which the ground never freezes, and Scenario B, where the ground freezes seasonally. For each scenario, the simulated ERT data was temperature corrected using the methodology outlined in the previous section (after Hayley et al., 2010). With these two scenarios, we aim to evaluate how well the temperature correction procedure performs under unfrozen and partially frozen conditions.

Resistivity model parameters
Assuming the rock physics model described in Equation 6 and the SFCC in Equation 7, we require two scenario-specific, time varying parameters to create a 3D bulk resistivity model: pore fluid salinity and subsurface temperature. Additionally, we require fitting parameters for Equations 6 and 7 that remain constant throughout the scenario. These resistivity model parameters are described in detail below.

Fluid salinity
Frack ponds contain highly saline water, with TDS concentrations that can exceed 300 g L −1 (Fakhru'l-Razi et al., 2009). Moderately saline flowback fluid has TDS concentrations of about 150 g L −1 (Gregory et al., 2011;Haluszczak et al., 2013), and we chose to use this representative value in this experiment. Background fluid TDS was assigned a value of 0.5 g L −1 , which would be expected of fresh water (Government of Canada, 2009).

Petrophysical relationships
In this experiment, we modeled the subsurface as a homogeneous, fully saturated fine-grained sand, so the estimated SFCC and best-fit Archie parameters for a fully saturated finegrained sandy material from the laboratory experiment carried out by Herring et al. (2019) were used in this study. Although simplified, this earth model is a reasonable starting point for a synthetic modeling experiment. The variables used to calculate bulk resistivity using Equation 6 are listed in Temp. compensation factor, ,˚C −1 0.0215 F I G U R E 2 Bulk resistivity as a function of temperature for modeled background conditions (i.e., where the total dissolved solids [TDS] concentration is equal to 0.5 g L −1 , as described in a previous laboratory experiment [Herring et al., 2019]). The orange inset shows a magnified version of the curve at temperature >0˚C Table 1, and the predicted bulk resistivity-temperature curve for background conditions is shown in Figure 2.

Subsurface temperature
We generated subsurface temperature profiles using a onedimensional (1D) conductive heat transport model: where T is temperature (˚C), t is time (s), z is depth (m), and α is thermal diffusivity (m 2 s −1 ).

Scenario A: Unfrozen ground
To create the subsurface temperature models for Scenario A, we used a sinusoidal surface temperature function that ranged from 0.1 to 30.0˚C (note that this surface temperature curve has the same amplitude as that in Scenario B). The thermal diffusivity of the soil was set to 1 × 10 −7 m 2 s −1 (Hinkel, 1997). The surface temperature and thermal diffusivity were used to generate 1D subsurface temperatures using Equation 18 (Figure 3a). The temperature modeling mesh used fine vertical discretization (0.01 m) to produce high-resolution subsurface temperature profiles, which were then interpolated onto the coarser mesh used for the resistivity model. Combining the resulting subsurface temperature, fluid salinity, and necessary fitting parameters in Equation 6 yielded the bulk resistivity model shown in Figure 3b. Because the subsurface never froze in Scenario A, resistivity variations due to seasonal temperature were smooth and low magnitude, similar to other studies where temperature corrections have been used successfully.

Scenario B: Seasonally frozen ground
In Scenario B, we created an earth model that was representative of conditions in northern British Columbia, where most Canadian fracking operations occur (Rivard et al., 2014). Monthly air temperature data from Fort St. John (Government of Canada, 2019) were fit with a sinusoid that ranged from −12.9 to 17.0˚C (Figure 4a). The temperature model for Scenario B also included a snow layer, where snow depths were interpolated from average monthly values (Government of Canada, 2019) with a spline function (Figure 4b). Thermal diffusivities of the snow layer, frozen soil, and unfrozen soil were chosen to be 1 × 10 −7 , 1 × 10 −6 , and 1 × 10 −7 m 2 s −1 , respectively, reflecting typical literature values (Hinkel, 1997;Liu & Si, 2008;Oldroyd et al., 2013). The resulting temperature models are shown in Figure 5a, and the bulk resistivity models are shown in Figure 5b. Note that the effects of latent heat were not considered in this temperature modeling. Consequently, the frozen layer is colder and more resistive in this synthetic scenario than it would be under real field conditions. For the two scenarios presented here, the phase and amplitude of the input air temperatures, the fluid salinity, and petrophysical relationships are the same, but because the air temperatures in Scenario B are 13˚C lower, the resulting resistivity models are very different. In Scenario A, the temperature related resistivity fluctuations are smooth and gradual, and in Scenario B, the seasonal freezing causes a large, sharp increase in resistivity in the frozen layer.

Vadose Zone Journal
F I G U R E 3 (a) Subsurface temperature profiles and (b) bulk resistivity models for Scenario A, where the no freezing occurs. Cross-sections are taken at the midpoint of the electrical resistivity tomography (ERT) survey and the red line shows the plane of the ERT survey. Note that this plot shows every other time step (i.e., every 2 mo). t is time, T is temperature, z is depth, y is lateral distance, and ρ is resistivity

Simulating ERT data
Forward modeling was used to simulate the ERT data that would be observed at monthly intervals for the duration of the monitoring experiment. We chose to model an ERT survey that used 48 electrodes and 5-m electrode spacing, so that the depth of investigation would be appropriate to image the plume. An ERT survey's ability to resolve features in the subsurface also depends on which electrode configurations are used. Two of the most commonly used field configurations, the Wenner array and the dipole-dipole array, were used here.
The Wenner array has a larger depth of penetration and is less sensitive to noise, whereas the dipole-dipole array has improved lateral resolution . We generated data using these two survey types to produce a more comprehensive dataset that takes advantage of the strengths of both array types. There were 315 measurements in the Wenner dataset and 990 in the dipole-dipole dataset, for a total of 1,305 data points. Because the data were generated using a discretized approach, choosing an appropriate 3D mesh for computa-tional modeling was an important part of minimizing numerical errors. In this case the minimum cell size was 0.17 × 0.25 × 0.25 m at the survey line, and cell sizes expanded with depth and lateral distance from the survey. When forward modeling a homogenous earth discretized using this mesh, apparent resistivities had <1% error. In synthetic modeling studies, it is standard practice to introduce normally distributed noise to the data to simulate field conditions, and prior studies have used between 1 and 5% Gaussian noise (Bastani et al., 2012;Miller et al., 2008;Morelli & Labrecque, 1996). For this work, where we are seeking to understand what we can achieve in a best-case scenario, we chose to add 1% noise. Noisy data noise (Ω) were created as follows: where no_noise (Ω) is the noise-free dataset and is a vector of normally distributed random noise. All ERT forward modeling and inversion was done using SimPEG, an open-source Python-based geophysical modeling framework (Cockett et al., 2015).

ERT data processing
The ERT data were inverted in 2D using the nonlinear least-squares approach described above. The regularization requires a few important inputs: the value of the regularization parameter β, the coefficients λ , λ , and λ that weight the model smallness and lateral and vertical smoothness constraints, and the reference model ref .
The simplest way of choosing β in practice is to fix it at a value that produces the smoothest model that still fits the observed data to within the specified noise level (Pidlisecky et al., 2007). We tested a range of β values between 0.5 and 500 and chose to use a fixed value of 5 for all inversions. For the initial inversions, λ and λ were set to 1 and λ was set to 10 −6 . ref was set to a homogenous model with the average apparent resistivity of each dataset, which is standard practice in the absence of prior structural information (Arato et al., 2015).

Vadose Zone Journal
F I G U R E 5 (a) Subsurface temperature profiles and (b) bulk resistivity models for Scenario B, where the ground freezes seasonally. Cross-sections are taken at the midpoint of the electrical resistivity tomography (ERT) survey and the red line shows the plane of the ERT survey. Note that this plot shows every other time step (i.e., every 2 mo) . t is time, T is temperature, z is depth, y is lateral distance, and ρ is resistivity The inversion process was ended when one of three convergence criteria was met: (a) the desired data misfit was achieved, (b) the objective function gradient decreased below a specified threshold, or (c) the maximum number of iterations was reached. Following Li and Oldenburg (2000), the target data misfit was set to Φ d = 1/2N, where N is the number of ERT data points (note that the 1/2 coefficient follows from the objective function formulation). The tolerance on the gradient of the objective function was chosen to be 0.1, and the maximum number of iterations was set to six. Although not all the iterations were able to fit the observed data to the specified noise levels, this is unsurprising as the inversion is attempting to fit the data with a 2D resistivity model, which is inherently unable to accurately reproduce the out-of-plane effects of the original 3D earth model. Inputting a higher maximum number of iterations increased computation time but produced negligible improvements in the data misfit, so an upper bound of six iterations is reasonable for this experiment.
Temperature corrections were applied to the data using the Hayley et al. (2010) method as described above. The temperature-resistivity relationship from Equation 16 was used to correct the estimated model est to TC est . We chose to correct the estimated models to a reference temperature of 15˚C in Scenario A and 2˚C in Scenario B. These temperatures are representative of the conditions at depth, beyond the influence of seasonal variations. This strategy minimizes the magnitude of the overall temperature correction, which is also helpful in field conditions where there is likely error in the parameterization of the resistivity-temperature relationship. In this experiment, we assume the petrophysical relationships and associated parameters are known exactly, as error analysis in this context is beyond the scope of this paper.

Scenario A: Plume detection in unfrozen ground
To assess the success of the temperature correction we chose to plot the apparent resistivities for midpoints at the center of the survey line. We did this by averaging the two to three apparent resistivity data points with midpoints nearest to = 0 for each pseudodepth for the Wenner and dipole-dipole arrays. Plotting these central apparent resistivities over time gives a holistic view of how seasonal temperature variations affect the data and how well the methodology outlined here can remove the effects of temperature. Figure 6 shows the apparent resistivity data for Scenario A. In unfrozen ground, the influence of near-surface temperature variability on measured apparent resistivities is relatively small and only affects the shallowest data points (Figure 6b). The temperature correction effectively removes the effect of temperature variations (Figure 6c), and the temperature corrected apparent resistivities contain less than 1.7% error (Figure 6d). The highest error is present in the shallowest data points, where the temperature correction underestimates the true apparent resistivity. At larger offsets the error in the temperature corrected data is minimal.
To make more quantitative assessments, we selected representative "shallower" (pseudodepth = 3.3) and "deeper" (pseudodepth = 10) data points from the dipole-dipole dataset and plotted them through time (Figure 7). The shallower pseudodepth represents a region where resistivity variation is primarily due to temperature, and the deeper pseudodepth represents a region where data are most influenced by the plume. Over the entire monitoring period, the mean absolute errors of these representative shallower and deeper data points are 0.79 and 0.10%, respectively (Table 2). In this scenario where the ground remains unfrozen, the temperature correction effectively reduces the deviations in apparent resistivity values caused by temperature variability.

Scenario B: Plume detection in seasonally frozen ground
In Scenario B where the shallow subsurface freezes seasonally, there is a large increase in apparent resistivities for small pseudodepths when the ground freezes, as a result of the Vadose Zone Journal F I G U R E 6 Scenario A apparent resistivities (ρ a ) at central midpoints over the 1-yr monitoring period for (a) plume + noise, (b) plume + noise + temperature, and (c) temperature corrected data. Panel d shows the normalized error between Panels a and c frozen surface layer (Figure 8b). The influence of temperature is also seen in larger offset data points, particularly for data collected in May (t = 334 d). This time step appears anomalous relative to adjacent time steps because it is the only sampled time step where there is a thawed layer (1-m thickness) overlaying the frozen layer (1-m thickness). Since electrical current tends to be concentrated in conductive regions of the subsurface and resistive features tend to have a shielding effect , even very large offsets have high apparent resistivities due to the conductive-resistive-conductive layering during this time step. This also demonstrates that while pseudodepth is often used as a proxy for depth, the two are not equivalent because sensitivity is affected by subsurface resistivity.
The temperature correction in partially frozen ground (Figure 8c) clearly does not perform as well as it did in Scenario A (Figure 6c). For most data points collected under partially frozen conditions, the temperature correction underestimates the true apparent resistivity by up to 79% (Figure 8d). The temperature corrected data contains significant error when the ground is partially frozen, and error in the shallower data is roughly proportional to the thickness of the frozen layer (Figure 9). Mean absolute error over the entire monitoring period in the shallower data is 29 and 20% in the deeper data ( Table 2). The reasons for the high error in the temperature corrected data, and attempts to reduce error through standard data processing and survey design strategies, are discussed below.

F I G U R E 7
Results from Scenario A showing the temporal variation of (a) temperature at various depths z, (b) distance between the plume's center of mass and the plane of the survey, and apparent resistivity (ρ a ) data for representative (c) shallower and (d) deeper data points. Shallower and deeper data points correspond to rows of the dipole-dipole plots in Figure 6 at pseudodepths of 3.3 and 10, respectively

ANALYSIS
It is clear that the ERT data processing strategies that work well to remove the effects of temperature in unfrozen ground are not readily transferrable to partially frozen ground. However, there are a few options in how ERT data is acquired and inverted that may improve how well we can resolve (and remove the influence of) a highly resistive surface layer. In the section below, we test different regularization strategies and electrode geometries to see whether an improvement in the inversion results translates to an improvement in the temperature-corrected data. For reference, the inverted models for = 304 d (see Figure 5b), when the plume's center of mass is 6.7 m from the plane of the ERT survey and the frozen layer is 2 m thick, are shown in Figure 10. We present a summary of error metrics for the temperature corrected representative shallower and deeper data points in Figure 11 and Table 2.

Higher lateral smoothness constraint
Regularization parameters, like the lateral smoothness constraint α , can be used to incorporate a priori knowledge of subsurface structure into the inversion process. In frozen ground, we expect a highly resistive surface layer that is F I G U R E 8 Scenario B apparent resistivities (ρ a ) at central midpoints over the 1-yr monitoring period for (a) plume + noise, (b) plume + noise + temperature, and (c) temperature corrected data. Panel d shows the normalized error between Panels a and c. Note that the color scale for Panels a-c has been clipped so as to not saturate the images laterally continuous, so a higher λ is appropriate to promote smoothness in the x direction. We tested a range of λ values on a representative time step where the top layer was frozen and found that λ = 500 produced a model with a smoother top layer without increasing the data misfit or smearing the plume appreciably when compared with the original results using λ = 1. We then followed the same ERT data processing steps described above using λ = 500. Although the estimated models looked smoother and more geologically reasonable (Figure 10c), the results of the temperature correction ( Figure 11) showed no noticeable improvement, with mean absolute errors equal to 28 and 16% for the shallower and deeper data, respectively (Table 2).

Disconnect inversion
Tikhonov regularization typically promotes smoothness between adjacent cells, but by relaxing this constraint between specific cells, it is possible to allow for larger changes across known boundary locations (Elwaseif & Slater, 2012). In this case, we multiplied λ by an identity matrix with zeros at frozen-unfrozen interface locations. We also tried using higher values of λ in order to enforce vertical smoothness within units. Using a value of λ = 5 within units and λ = 0 in boundary locations produced the most geologically reasonable estimated model (Figure 10d), but the frozen layer was still not well resolved. This could potentially be due to F I G U R E 9 Results from Scenario B showing the temporal variation of (a) temperature at various depths z, (b) frozen layer thickness and distance between the plume's center of mass and the plane of the survey, and apparent resistivity (ρ a ) data for representative (c) shallower and (d) deeper data points. Shallower and deeper data points correspond to rows of the dipole-dipole plots in Figure 8 at pseudodepths of 3.3 and 10, respectively large electrode spacing relative to the thickness of the frozen surface layer. Ultimately, the temperature correction was not improved from Scenario B (see Figure 11 and Table 2).

L1 model misfit norm
An L1 measure of model misfit can be used to reduce smoothing of features and damping of amplitudes in the inversion process. We tested an IRLS approach where the maximum number of iterations was set to 40 and the maximum number of IRLS iterations was set to 20. Using a higher number of iterations is common practice, as this method requires more iterations to converge (Farquharson & Oldenburg, 1998). This method produced some artifacts (Figure 10e) as it attempted to more closely fit the data noise, particularly for the dipoledipole array that is more sensitive to noise. The results of the temperature correction in this case did not improve significantly ( Figure 11, Table 2).

Burying electrodes beneath frozen layer
One option to minimize the effect of temperature on ERT data is to bury electrodes below the maximum depth of freezing. We simulated data with electrodes buried below the frost zone at a depth of 2.5 m. The temperature correction using buried electrodes performed quite poorly. Temperature corrected apparent resistivities had a larger error when the ground was unfrozen and were overestimated by up to 30% in the shallow data points. Temperature corrected apparent resistivities were also underestimated when freezing occurred. The source of these errors was apparent in the inverted resistivity model (Figure 10f). The data collected with buried electrodes was sensitive to the ground all around the electrodes and lacked the directional information that could localize temperature related features above the electrodes. The estimated resistivity models therefore produced worse estimates of subsurface F I G U R E 1 0 Estimated resistivity models for time (t) = 304 d (see Figure 5b) for (a) Scenario A (no freezing), (b) Scenario B (seasonal freezing), (c) Scenario B inverted with a higher lateral smoothness constraint (λ = 500), (d) Scenario B inverted with a disconnect inversion, (e) Scenario B inverted with a L1 model misfit norm, (f) Scenario B with data collected using buried electrodes, and (g) Scenario B with data collected using a smaller electrode spacing (1 m). Note that the color scale is different in Panel a, and lateral and vertical limits are different in Panel g. x denotes horizontal distance along the survey line, z is depth, and ρ is resistivity in Ω m resistivity, causing the temperature correction to perform poorly ( Figure 11, Table 2).

Smaller electrode spacing
As noted earlier, standard practices for resolving sharp boundaries did not significantly improve the temperature-corrected data. We hypothesized that the relatively large electrode spacing (5 m) compared with the thickness of the frozen layer (maximum ∼2 m) led to poor resolution of the frozen layer, regardless of the data processing strategy used. A reasonable next step was to test an alternative ERT survey design that used 120 electrodes at 1-m spacing. We also simulated a highdensity dataset, where the maximum a spacing was set to 31 F I G U R E 1 1 Comparison of the error in temperature corrected data for all test cases. Scenario A and Scenario B use a surface electrical resistivity tomography (ERT) array with 48 electrodes at 5-m electrode spacing, λ = 1, and an L2 norm. The shallower data points shown here use quadripoles with small pseudodepths and midpoints at the center of the array, and are most sensitive to temperature effects (Panel a). Deeper data use higher pseudodepths and are more sensitive to the plume (Panel b) and the maximum n value was set to 10. This produced a dense dataset with 13,069 data points at each time step. The downside of this setup is that data acquisition and processing time is much higher than with the previous array. Since the pseudodepths were different for this experiment, we picked values of 0.67 and 8.67 as representative shallower and deeper data points, respectively. Although the frozen layer is better resolved in the inverted model, the depth of investigation is lower and the resolution of the plume is poor (Figure 10g). Interestingly, the shallow temperature-corrected data had a much higher misfit than previous data processing and acquisition strategies ( Figure 11); on average, about 805% ( Table 2). The reasons for this misfit will be explored in the section below.

DISCUSSION
Although marginal improvements can be made in the temperature correction by adjusting how ERT data were processed, the effects of temperature were not corrected satisfactorily under frozen conditions, even in this idealized scenario. In the section above, we showed that changing how data are collected, either by burying the electrodes or by using smaller electrode spacing, can actually increase errors in the temperature-corrected data.
F I G U R E 1 2 Schematic for two-layered medium. C1 and P1 are current and potential electrodes. T is temperature (˚C), z is depth (m), and ρ is resistivity (Ω m) With these results in mind, we hypothesized that the temperature correction fails under partially frozen conditions either (a) because the magnitude of the temperature correction Δ is too large and the first-order Taylor approximation is not appropriate, or (b) because the resistivity of the frozen layer is underestimated in the estimated model, leading to errors in the temperature correction. To test this, we performed a sensitivity analysis to examine how error in the temperature corrected data is affected by the magnitude of the temperature correction as well as error in the estimated model. This sensitivity analysis was carried out using a simple two-layered medium where the homogeneous resistivities of top and bottom layers ρ 1 and ρ 2 (Ω m) depend on their temperatures 1 and 2 (˚C) (Figure 12). We then forward modeled ERT data and applied a temperature correction. By giving the frozen surface layer F I G U R E 1 3 (a) Maximum absolute and (b) mean absolute errors in temperature corrected electrical resistivity tomography (ERT) data generated using a two-layered analytical solution. Error in the estimated model was introduced by underestimating the resistivity of the top layer (ρ 1 ). Inlays show data errors for when the temperature of the top layer 1 ≥ 0˚C. Note that the y axis differs between plots. The ERT data were generated using an array with 5-m electrode spacing (like the one used for Scenarios A and B) and an array with 1-m spacing (as described in Section 5) a range of resistivity values, we tested how the results were affected by the magnitude of the temperature correction. By introducing a range of error into the estimated ρ 1 , we tested how sensitive the temperature correction was to errors in the inverted model.
To approximate Scenario B, the depth to the interface (m) was set to 2 m (the maximum depth of freezing), below which the temperature was a uniform 2˚C. To calculate the resistivity of each layer, we used the same petrophysical relationships as before. An analytical solution exists for the potential field (V) due to a point source at the surface for this simple layered earth (Telford et al., 1990): where (A) is the current, (m) is the distance between the point source and the measurement location, and (unitless) is a reflection coefficient: We used the analytical solution here because it minimized computation time and avoided the numerical errors associated with numerical modeling.
For this sensitivity testing, the surface layer was assigned a range of temperatures between −20 and 20˚C. Error was introduced into the "estimated" model by underestimating the resistivity of the top layer, ρ 1 . The temperature-corrected data were calculated the same way as before to a reference temperature of 2˚C. Figure 13 shows how the mean and maximum errors in temperature corrected data change as a function of surface layer temperature, electrode spacing, and error in ρ 1 . The results from this sensitivity analysis show that the temperaturecorrected data are much more sensitive to errors in the estimated model under partially frozen conditions than it is under unfrozen conditions. When the top layer is frozen, even a 2% underestimation of ρ 1 causes a maximum of 9-29% error in the temperature corrected data for an ERT survey with 5-m electrode spacing. It is also worth noting that the data points that experience the maximum errors are those with the smallest offsets. These findings illustrate that the data temperature correction of Hayley et al. (2010) is indeed sensitive to error in the inverted resistivity model, and although this sensitivity is minimal above 0˚C, it can greatly affect temperature corrections at subzero temperatures.
Another key finding is that the survey with 1-m electrode spacings is much more sensitive to errors in the estimated model, compared with the data collected using 5-m electrode spacings. This helps to explain why using smaller electrode spacing decreased the accuracy of the temperature corrected data collected at short offsets. Although a more refined ERT survey can improve resolution of the frozen near-surface layer, the increased sensitivity caused a larger error in shallow temperature corrected data. This is an important caveat to consider when designing field surveys.
Interestingly, when there is no error in the estimated model, the temperature correction also works perfectly, regardless of the magnitude of Δρ. This analysis demonstrates that the limitations in using this method to correct for the effects of a seasonally frozen surface layer are mostly due to errors in the estimated model, rather than issues with the first-order approximation in the temperature correction method.

CONCLUSIONS
In this study, we demonstrated that even under ideal conditions, we cannot satisfactorily remove the influence of a seasonally frozen surface layer from the data collected in a time-lapse ERT monitoring experiment. The ERT data collection and processing strategies that are often used in frozen ground, although usually producing more geologically reasonable inverted models, did not appreciably improve the outcome of the temperature correction and, in some cases, caused the error in the temperature corrected data to increase. A sensitivity analysis showed that errors in the temperaturecorrected data were primarily due to error in the estimated resistivity model, and that this sensitivity to model error was much higher in partially frozen ground than in unfrozen ground. Ultimately, successfully removing the effects of temperature from ERT data collected in frozen ground requires a more accurate estimate of subsurface resistivity than is achieved with current best practice data collection and processing methods. Whether or not the target of the ERT survey can be resolved also depends on the magnitude of the signal of interest relative to the signal from temperature variation, and this relationship could be explored in future work. In order to image and remove the effects of a frozen layer, future studies may need to reframe the ERT inversion process entirely to avoid the limitations that are inherent to a pixel parameterization with a large number of model parameters. The findings from this study have important implications for any situation where the practitioner is trying to compensate for the effects of a large near-surface resistivity perturbation, whether due to freezing, drying, or wetting, or a changing water layer in marine ERT surveys. In these situations, it is important to acknowledge how uncertainty in the estimated model propagates into the corrected data, and ultimately how well surface perturbations can be removed to isolate the signal of interest.

A C K N O W L E D G M E N T S
This research was funded by a Mitacs Accelerate Ph.D. fellowship in partnership with Matrix Solutions, the University of Calgary, and the Natural Sciences and Engineering Research Council of Canada. Thanks to the associate editor Dr. Andrew Binley, Dr. Florian Wagner, and an anonymous reviewer for their feedback that greatly improved the quality of the manuscript.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.