Sensitivity of airborne geophysical data to sublacustrine and near-surface permafrost thaw

A coupled hydrogeophysical forward and inverse modeling approach is developed to illustrate the ability of frequency-domain airborne electromagnetic (AEM) data to characterize subsurface physical properties associated with sublacustrine permafrost thaw during lake-talik formation. Numerical modeling scenarios are evaluated that consider non-isothermal hydrologic responses to variable forcing from different lake depths and for different hydrologic gradients. A novel physical property relationship connects the dynamic distribution of electrical resistivity to ice saturation and temperature outputs from the SUTRA groundwater simulator with freeze–thaw physics. The influence of lithology on electrical resistivity is controlled by a surface conduction term in the physical property relationship. Resistivity models, which reflect changes in subsurface conditions, are used as inputs to simulate AEM data in order to explore the sensitivity of geophysical observations to permafrost thaw. Simulations of sublacustrine talik formation over a 1000-year period are modeled after conditions found in the Yukon Flats, Alaska. Synthetic AEM data are analyzed with a Bayesian Markov chain Monte Carlo algorithm that quantifies geophysical parameter uncertainty and resolution. Major lithological and permafrost features are well resolved by AEM data in the examples considered. The subtle geometry of partial ice saturation beneath lakes during talik formation cannot be resolved using AEM data, but the gross characteristics of sub-lake resistivity models reflect bulk changes in ice content and can identify the presence of a talik. A final synthetic example compares AEM and ground-based electromagnetic responses for their ability to resolve shallow permafrost and thaw features in the upper 1–2 m below ground outside the lake margin.


Introduction
Permafrost thaw can have important consequences for the distribution of surface water (Roach et al., 2011;Rover et al., 2012), stream discharge and chemistry (O'Donnell et al., 2012;Petrone et al., 2007;Striegl et al., 2005;Walvoord and Striegl, 2007), and exchange between groundwater and surface water systems (Bense et al., 2009;Callegary et al., 2013;Walvoord et al., 2012).Likewise, hydrologic changes that alter the thermal forcing supplied by surface water or groundwater systems can modify the distribution of permafrost, illustrating the strong feedbacks between permafrost and hydrology.In addition to hydrologic processes, permafrost is affected by climate warming in Arctic and sub-Arctic regions (Hinzman et al., 2005;Jorgenson et al., 2001) as well as disturbance by fire (Yoshikawa et al., 2002).Climate feedbacks associated with permafrost thaw include changes in the amount of organic carbon stored in soils that is vulnerable to decomposition (Koven et al., 2011;O'Donnell et al., 2011) and subsequent methane and carbon dioxide released from soils by the degradation of organic material previously sequestered in frozen ground (Anthony et al., 2012).Permafrost thaw also has significant implications for land management and infrastructure, including the potential to damage buildings, roadways, or pipelines due to ground settling, and thermal erosion that can alter coastlines and landscape stability (Larsen et al., 2008;Nelson et al., 2002).
Several investigations have shown the significance of climate and advective heat transport in controlling the distribution of permafrost in hydrologic systems (Bense et al., 2009;Rowland et al., 2011;Wellman et al., 2013).These results yield important insight into the mechanistic behavior of coupled thermal-hydrologic systems and are a means for predicting the impact on permafrost from a wide range of climate and hydrologic conditions.However, few techniques are capable of assessing the distribution of permafrost, and most approaches only capture a single snapshot in time.
Satellite remote-sensing techniques have proven useful in detecting the distribution and changes in shallow permafrost, vegetation, and active layer thickness over large areas (Liu et al., 2012;Panda et al., 2010;Pastick et al., 2014) but are only sensitive to very near-surface properties.Borehole cores and downhole temperature or geophysical logs provide direct information about permafrost and geologic structures but tend to be sparsely located and are not always feasible in remote areas.Geophysical methods are necessary for investigating subsurface physical properties over large and/or remote areas.Recent examples of geophysical surveys aimed at characterizing permafrost in Alaska include an airborne electromagnetic (AEM) survey used to delineate geologic and permafrost distributions in an area of discontinuous permafrost (Minsley et al., 2012), ground-based electrical measurements used to assess shallow permafrost aggradation near recently receded lakes (Briggs et al., 2014), electrical and electromagnetic surveys used to characterize shallow active layer thickness and subsurface salinity (Hubbard et al., 2013), and surface nuclear magnetic resonance soundings used to infer the thickness of unfrozen sediments beneath lakes (Parsekian et al., 2013).A challenge with geophysical methods, however, is that geophysical properties (e.g., electrical resistivity) are only indirectly sensitive to physical properties of interest (e.g., lithology, water content, thermal state).In addition, various physical properties can produce similar electrical resistivity values.Therefore, it is critically important to understand the relationship between geophysical properties and the ultimate physical properties and processes of interest (Minsley et al., 2011;Rinaldi et al., 2011).
The non-isothermal hydrologic simulations of Wellman et al. (2013) predict the evolution of lake taliks (unfrozen sublacustrine areas in permafrost regions) in a two-dimensional axis-symmetric model under different environmental scenarios (e.g., lake size, climate, groundwater flow regime).Here, we investigate the ability of geophysical measurements to recover information about the underlying spatial distribution of permafrost and hydrologic properties.This is accomplished in three steps: (1) development of a physical property relation that connects permafrost and hydrologic properties to geophysical properties, (2) generation of synthetic geophysical data that would be expected for various permafrost hydro-logic conditions that occur during simulated lake-talik formation, and (3) inversion of the synthetic geophysical data using realistic levels of noise to investigate the ability to resolve specific physical features of interest.Our focus is on electromagnetic geophysical methods as these types of data have previously been acquired near Twelvemile Lake in the Yukon Flats, Alaska (Ball et al., 2011;Minsley et al., 2012); this lake is also the basis for the lake simulations discussed by Wellman et al. (2013).Wellman et al. (2013) describe numerical simulations of lake-talik formation in watersheds modeled after those found in the lake-rich Yukon Flats of interior Alaska.Modeling experiments used the SUTRA groundwater modeling code (Voss and Provost, 2002) enhanced with capabilities to simulate freeze-thaw processes (McKenzie and Voss, 2013).The phase change between ice and liquid water occurs over a specified temperature range and accounts for latent heat of fusion as well as changes in thermal conductivity and heat capacity for ice-water mixtures.Ice content also changes the effective permeability, thereby altering subsurface flowpaths and enforcing a strong coupling between hydraulic and thermal processes.The modeling domain, which is adapted for this study, is axis-symmetric with a central lake and upwards-sloping ground surface that rises from an elevation of 500 m at r = 0 to 520 m at the outer extent of the model, r = 1800 m (Fig. 1).The model uses a layered geology consistent with the Yukon Flats (Minsley et al., 2012;Williams, 1962), with defined hydrologic and geophysical parameters for each layer summarized in Table 1.Initial permafrost conditions prior to lake formation were established by running the model to steady state under hydrostatic conditions with a constant temperature of −2.25 • C applied to the land surface, which produces a laterally continuous permafrost layer extending to a depth of about 90 m.

Coupled thermal-hydrologic simulations
Subsequent hydrologic simulations assume fully saturated conditions and are performed over a 1000-year period under 36 different scenarios of climate (warmer than, colder than, and similar-to-present conditions), hydrologic gradient (hydrostatic, gaining, and losing lake conditions), and lake depth/extent (3, 6, 9, and 12 m deep lakes that intersect the ground surface at increasing distance, as shown in Fig. 1).Complete details and results of the hydrologic simulations can be found in Wellman et al. (2013).At each simulation time step, the SUTRA model outputs temperature, pressure, and ice saturation.Conversion of these hydrologic variables to electrical resistivity -the geophysical property needed to simulate electromagnetic data considered here -is described below.

A physical property relationship
Electrical resistivity is the primary geophysical property of interest for the electromagnetic geophysical methods used in this study.It is well established that resistivity is sensitive to basic physical properties such as unfrozen water content, soil or rock texture, and salinity (Palacky, 1987).Here, we build on earlier efforts to simulate the electrical properties of ice-saturated media (Hauck et al., 2011) by using a modified form of Archie's Law (Archie, 1942) that also incorporates surface conduction effects (Revil, 2012) to predict the dynamic electrical resistivity structure for the evolving state of temperature and ice saturation (S i ) in the talik simulations.Bulk electrical conductivity is described by Revil (2012) as where σ is the bulk electrical conductivity [S m −1 ]; S w is the fractional water saturation [-] in the pore space, where S w = 1−S i ; σ f is the conductivity of the saturating pore fluid [S m −1 ]; m is the Archie cementation exponent [-]; n is the Archie saturation exponent [-]; F is the formation factor [-], where F = φ −m and φ is the matrix porosity [-]; and σ s is the conductivity [S m −1 ] associated with grain surfaces.The Archie exponents m and n are known to vary as a function of pore geometry; here, we use m = n = 1.5, which is appropriate for unconsolidated sediments (Sen et al., 1981).Simulation results are presented as electrical resistivity [ m], which is the inverse of the conductivity, i.e., ρ = 1/σ .The first term in Eq. (1) describes electrical conduction within the pore fluid, where fluid conductivity is defined as (2) The summation in Eq. ( 2) is over all dissolved ionic species (Na + and Cl − are assumed to be the primary constituents in this study), where F c is Faraday's constant [C mol −1 ] and C i , β i , and z i are the concentration [mol L −1 ], ionic mobility [m 2 V −1 s −1 ], and valence of the ith species, respectively.1.5 Water salinity (C) [ppm] 250 200, 10, 500 Salinity exponent (a) [-] 0.8 Surface conduction effects, described by the second term in Eq. ( 1), are related to the chemistry at the pore-water interface, and can be important in fresh water (low conductivity) systems at low porosity (high ice saturation).Additionally, the surface conduction term provides a means for describing the conductivity behavior for different lithologies, as will be described below.The surface conductivity is given by where β s is the cation mobility [m 2 V −1 s −1 ] for counterions in the electrical double layer at the grain-water interface (Revil et al., 1998) and Q v is the excess electrical charge density [C m −3 ] in the pore volume, and where ρ g is the mass density of the grains [kg m −3 ] and χ is the cation exchange capacity [C kg −1 ].Changes in χ, representative of bulk differences in clay mineral content, are used to differentiate the electrical signatures of the lithologic units in this study (Table 1).
The temperature, T [C], dependence of ionic mobility affects both the fluid conductivity (Eq.2) and surface conductivity (Eq.3), where mobility is approximated as a linear function of temperature (Keller and Frischknecht, 1966;Sen and Goode, 1992) as (5) Finally, we consider the effect of increasing ice saturation on salinity.Because salts are generally excluded as freezing occurs, salinity of the remaining unfrozen pore water is expected to increase with increasing ice content (Marion, www.the-cryosphere.net/9/781/2015/The Cryosphere, 9, 781-794, 2015 1995), leading to a corresponding increase in fluid conductivity according to Eq. (2).To describe this dependence of salinity on ice saturation, C(S i ), we use the expression where α ∼ 0.8 accounts for loss of solute from the pore space due to diffusion or other transport processes and S i = 1−S w .
Information about the different lithologic units described by Wellman et al. (2013) that are also summarized in Table 1 are used to define static model properties such as porosity, grain mass density, cation exchange capacity, and Archie's exponents.Dynamic outputs from the SUTRA simulations, including temperature and ice saturation, are combined with the static variables in Eqs. ( 1)-( 6) to predict the evolving electrical resistivity structure.

Geophysical forward simulations
Synthetic AEM data are simulated for each snapshot of predicted bulk resistivity values using nominal system parameters based on the Fugro RESOLVE 1 frequency-domain AEM system that was used in the Yukon Flats survey (Minsley et al., 2012).The RESOLVE system consists of five horizontal coplanar (HCP) transmitter-receiver coil pairs separated by approximately 7.9 m that operate at frequencies 0.378, 1. 843, 8.180, 40.650, and 128.510 kHz and one vertical coaxial coil pair with 9 m separation that operates at 3.260 kHz.Oscillating currents and associated magnetic fields created by the transmitter coils induce electrical currents in the subsurface that, in turn, generate secondary magnetic fields that are recorded by the receiver coils (Siemon, 2006;Ward and Hohmann, 1988).Data are reported as in-phase and quadrature components of the secondary field in parts per million (ppm) of the primary field, and responses as a function of frequency can be converted through mathematical inversion to estimates of electrical resistivity as a function of depth (e.g., Farquharson et al., 2003).Data are simulated at the nominal survey elevation of 30 m above ground surface using the one-dimensional modeling algorithm described in Minsley (2011), which follows the standard electromagnetic theory presented by Ward and Hohmann (1988).
The vertical profile of resistivity as a function of depth is extracted at each survey location and is used to simulate forward geophysical responses.There are 181 sounding locations for each axis-symmetric model, starting at the center of the lake (r = 0 m) and moving to the edge of the model domain (r = 1800 m) in 10 m increments.Each vertical resistivity profile extends to 200 m depth, which is well beyond the depth to which we expect to recover parameters in the geophysical inversion step.A center-weighted five-point filter with weights equal to [0.0625, 0.25, 0.375, 0.25, 0.0625] is used to average neighboring bulk resistivity values at each 1 Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government.
depth before modeling in order to partly account for the lateral sensitivity of AEM systems (Beamish, 2003).Forward simulations are repeated for each of the 50 simulation times between output of 0 and 1000 years from SUTRA, resulting in 9050 data locations per modeling scenario.
Synthetic ground-based electromagnetic data presented in Sect.3.3 are simulated using nominal system parameters based on the GEM-2 instrument (Huang and Won, 2003).The GEM-2 has a single HCP transmitter-receiver pair separated by 1.66 m, and data are simulated at six frequencies: 1.5, 3.5, 8.1, 19, 43, and 100 kHz.A system elevation of 1 m above ground is assumed, which is typical for this handcarried instrument.

Parameter estimation and uncertainty quantification
The inverse problem involves estimating subsurface resistivity values given the simulated forward responses and realistic assumptions about data errors.Geophysical inversion is inherently uncertain; there are many plausible resistivity models that are consistent with the measured data.In addition, the ability to resolve true resistivity values is limited both by the physics of the AEM method and the level of noise in the data.Here, we use a Bayesian Markov chain Monte Carlo (McMC) algorithm developed for frequency-domain electromagnetic data (Minsley, 2011) to explore the ability of simulated AEM data to recover the true distribution of subsurface resistivity values at 20-year intervals within the 1000-year lake-talik simulations.This McMC approach is an alternative to traditional inversion methods that find a single "optimal" model that minimizes a combined measure of data fit and model regularization (Aster et al., 2005).Although computationally more demanding, McMC methods allow for comprehensive model appraisal and uncertainty quantification.AEM-derived resistivity estimates for the simulations considered here will help guide interpretations of future field data sets, identifying the characteristics of relatively young versus established thaw under different hydrologic conditions.
The McMC algorithm provides comprehensive model assessment and uncertainty analysis and is useful in diagnosing the ability to resolve various features of interest.At every data location along the survey profile, an ensemble of 100 000 resistivity models is generated according to the Metropolis-Hastings algorithm (Hastings, 1970;Metropolis et al., 1953).According to Bayes' theorem, each model is assigned a posterior probability that is a measure of (1) its prior probability which, in this case, is used to penalize models with unrealistically large contrasts in resistivity over thin layers; and (2) its data likelihood, which is a measure of how well the predicted data for a given resistivity model match the observed data within data errors.A unique aspect of this algorithm is that it does not presuppose the number of layers needed to fit the observed data, which helps avoid biases due to assumptions about model parameterization.Instead, trans-dimensional sampling rules (Green, 1995;Sambridge et al., 2013) are used to allow the number of unknown layers to be one of the unknowns.That is, the unknown parameters for each model include the number of layers, layer interface depths, and resistivity values for each layer.
Numerous measures and statistics are generated from the ensemble of plausible resistivity models, such as the single most-probable model, the probability distribution of resistivity values at any depth, the probability distribution of where layer interfaces occur as a function of depth, and the probability distribution of the number of layers (model complexity) needed to fit the measured data.A detailed description of the McMC algorithm can be found in Minsley (2011).Finally, probability distributions of resistivity are combined with assumptions about the distribution of resistivity values for any lithology and/or ice content in order to make a probabilistic assessment of lithology or ice content, as illustrated below.

Electrical resistivity model development
Information about the different lithologic units described by Wellman et al. (2013) that are also summarized in Table 1 are used to define static model properties such as porosity, grain mass density, cation exchange capacity, and Archie's exponents.Dynamic outputs from the SUTRA simulations, including temperature and ice saturation, are combined with the static variables in Eqs. ( 1)-( 6) to predict the evolving electrical resistivity structure.The behavior of bulk resistivity as a function of ice saturation is shown in Fig. 2. Separate curves are shown for a range of χ (cation exchange capacity) values, which are the primary control in defining offset resistivity curves for different lithologies, where increasing χ is generally associated with more fine-grained material such as silt or clay.
For each of the 1000-year simulations, the static variables summarized in Table 1 are combined with the spatially and temporally variable state variables T and S i output by SU-TRA to predict the distribution of bulk resistivity at each time step using Eqs.( 1)-( 6).An example of SUTRA output variables for the 6 m deep gaining lake scenario at 240 years (the approximate sub-lake talik breakthrough time for that scenario) is shown in Fig. 3a and b, and the predicted resistivity for this simulation step is shown in Fig. 3c.The influence of different lithologic units is clearly manifested in the predicted resistivity values, whereas lithology is not overly evident in the SUTRA state variables.For a single unit, there is a clear difference in resistivity for frozen versus unfrozen conditions.Across different units, there is a contrast in resistivity when both units are frozen or unfrozen.Resistivity can therefore be a valuable indicator of both geologic and ice Bulk resistivity as a function of ice saturation using the physical properties defined for each of the lithologic units described in Table 1.content variability.However, there is also ambiguity in resistivity values as both unfrozen unit 2 and frozen unit 3 appear to have intermediate resistivity values of approximately 100-300 m (Fig. 3c) and cannot be characterized by their resistivity values alone.This ambiguity in resistivity can only be overcome by additional information such as borehole data or prior knowledge of geologic structure.Synthetic bulk resistivity values according to Eq. ( 1) are shown in Fig. 4 for the four different lake depths (3, 6, 9, and 12 m) at three different simulation times (100, 240, and 1000 years) output from the hydrostatic/current climate condition simulations.
Lithology and ice saturation are the primary factors that control simulated resistivity values (Fig. 2), though ice saturation is a function of temperature.The empirical relation between temperature and bulk resistivity is shown in Fig. 3d by cross-plotting values from Fig. 3b and c.Within each lithology resistivity is relatively constant above 0 • C, with a rapid increase in resistivity for temperatures below 0 • C.This result is very similar to the temperature-resistivity relationships illustrated by Hoekstra et al. (1975,   their ability to discern differences among very high resistivity values, as discussed later, this artifact does not significantly impact the results presented here.

Parameter estimation and uncertainty quantification
AEM data (not shown) are simulated for each of the electrical resistivity models (e.g., Fig. 4) using the methods described in Sect.2.3.The simulated data are then used to re- The overall pattern of different lithologic units and frozen/unfrozen regions is accurately depicted in Fig. 5, with two exceptions that will be discussed in greater detail: (1) the specific distribution of partial ice saturation beneath the lake before thaw has equilibrated (Fig. 5a and b) and (2) the shallow sand layer (unit 1) that is generally too thin to be resolved using AEM data.
A point-by-point comparison of true (Fig. 4f) versus predicted (Fig. 5c) resistivity values for the hydrostatic 6 m deep lake scenario at the simulation time 1000 years is shown in Fig. 6a.The cross-plot of true versus estimated resistivity values generally fall along the 1 : 1 line, providing a more quantitative indication of the ability to estimate the subsurface resistivity structure.Estimates of the true resistivity values for each lithology and freeze-thaw state (Fig. 6b) tend to be indistinct, appearing as a vertical range of possible values in Fig. 6a due to the inherent resolution limitations of inverse methods and parameter tradeoffs (Day-Lewis et al., 2005;Oldenborger and Routh, 2009).Although the greatest point density for both frozen and unfrozen silts (unit 3) falls along the 1 : 1 line, resistivity values for these components of the model are also often overestimated; this is likely due to uncertainties in the location of the interface between the silt and gravel units.This is in contrast with the systematic underestimation of frozen gravel resistivity values due to the inability to discriminate very high resistivity values using electromagnetic methods (Ward and Hohmann, 1988).Frozen sands (true log resistivity ∼ 2.8 in Fig. 6b) are also systematically overestimated in Fig. 6a; in this case, due to the inability to resolve this relatively thin resistive layer.
While useful, single "best" estimates of resistivity values at any location (Fig. 6) are not fully representative of the information contained in the AEM data and associated model uncertainty.From the McMC analysis of 100 000 models at each data location, estimates of the posterior probability density function (pdf) of resistivity are generated for each point in the model.Probability distributions are extracted from a depth of 15 m, within the gravel layer (unit 2), at one location where unfrozen conditions exist (r = 0 m) and a second location outside the lake extent (r = 750 m) where the ground remains frozen (Fig. 7a).Results from a depth of 50 m, within the silt layer (unit 3), are shown in Fig. 7b.With the exception of the frozen gravels, the resistivity of which tends to be underestimated, the peak of each pdf is a good estimate of the true resistivity value at that location.
Resistivity values are translated to estimates of ice saturation, which is displayed on the upper axis of each panel in Fig. 7, using the appropriate lithology curve from Fig. 2. Using the ice-saturation-transformed pdfs, quantitative inferences can be made about the probability of the presence or absence of permafrost.For example, the probability of ice content being less than 50 % is estimated by calculating the fractional area under each distribution for ice-content values www.the-cryosphere.net/9/781/2015/The Cryosphere, 9, 781-794, 2015 less than 0.5.Probability estimates of ice content less than 50 % and greater than 95 % for the four distributions shown in Fig. 7 are summarized in Table 2. High probabilities of ice content exceeding 95 % are associated with the r = 750 m location outside the lake extent, whereas high probability of ice content below the 50 % threshold are observed at r = 0 beneath the center of the lake.The pdfs for each lithology shown in Fig. 7 are end-member examples of frozen and unfrozen conditions.Within a given lithology, a smooth transition from the frozen-state pdf to the unfrozen-state pdf is observed as thaw occurs, with corresponding transitions in the calculated ice threshold probabilities.Further illustration of the spatial and temporal changes in resistivity pdfs are shown in Fig. 8.The resistivity pdf is displayed as a function of distance from the lake center at the same depths (15 and 50 m) shown in Fig. 7, corresponding to gravel (Fig. 8a, c, and e) and silt (Fig. 8b, d, and f) locations.High probabilities, i.e., the peaks in Fig. 7, correspond to dark-shaded areas in Fig. 8. Images are shown for three different time steps in the SUTRA simulation for the hydrostatic 6 m deep lake scenario: 100 years (Fig. 8a and b), 240 years (Fig. 8c and d), and 1000 years (Fig. 8e and f).Approximate ice-saturation values, translated from the ice versus resistivity relationships for each lithology shown in Fig. 2, are displayed on the right axis of each panel in Fig. 8, and true resistivity values are plotted as a dashed line.Observations from Fig. 8 include the following: 1. Outside the lake boundary, pdfs are significantly more sharply peaked (darker shading) for the gravel unit than the silt unit, suggesting better resolution of shallower resistivity values within the gravel layer.It should be noted, however, that this improved resolution does not imply improved model accuracy; in fact, the highest probability region slightly underestimates the true resistivity value.
2. Probability distributions for the silt layer track the true values but with greater uncertainty.
3. Inside the lake boundary, gravel resistivity values are not as well resolved compared with locations outside the lake boundary due to the loss of signal associated with the relatively conductive lake water.
4. Increasing trends in resistivity/ice saturation towards the outer extents of the lake are captured in the pdfs but are subtle.
5. Within the silt layer at early times before the talik is fully through-going (Fig. 8b and d), the AEM data are insensitive to which layer is present, hence the bi-modal resistivity distribution with peaks associated with characteristic silt and gravel values.This ambiguity disappears at later times when the low-resistivity unfrozen silt layer extends to the base of the unfrozen gravels, which is a more resolvable target (Fig. 8f).
A more detailed analysis of the changes in resistivity and ice saturation as a function of time, and for the differences between hydrostatic and gaining lake conditions, is presented in Fig. 9. Average values of resistivity/ice saturation within 100 m of the lake center are shown within the gravel layer at a depth of 15 m (Fig. 9a) and a depth of 50 m within the silt layer (Fig. 9b) at 20-year time intervals.Outputs are displayed for both 6 m deep hydrostatic and gaining lake scenarios.Thawing due to conduction occurs over the first ∼ 200 years within the gravel layer (Fig. 9a), with similar trends for both the hydrostatic and gaining lake conditions and no clear relationship to the talik formation times indicated as vertical lines.Conduction-dominated thaw is observed for the gravel layer in the gaining lake scenario because significant advection does not occur until after the thaw bulb has extended beneath the gravel layer.In the deeper silt layer (Fig. 9b), however, very different trends are observed for the hydrostatic and gaining lake conditions.Ice content decreases gradually as thawing occurs in the hydrostatic scenario, consistent with conduction-dominated thaw, reaching a minimum near the time of talik formation at 687 years (Wellman et al., 2013, Table 3).In contrast, there is a rapid loss in ice content in the gaining lake scenario resulting from the influence of advective heat transport as groundwater is able to move upwards through the evolving talik beneath the lake.This rapid loss in ice content begins after the gravel layer thaws and reaches a minimum near the 258-year time of talik formation for this scenario.These trends, captured by the AEM-derived resistivity models, are consistent with the plots of change in ice volume output from the SUTRA simulations reported by Wellman et al. (2013, Fig. 3).

Near-surface resolution
Finally, we focus on the upper sand layer (unit 1), which is generally too thin (2 m) and resistive (> 600 m) to be resolved using AEM data, though may be imaged using other ground-based electrical or electromagnetic geophysical methods.Seasonal thaw and surface runoff cause locally reduced resistivity values in the upper 1 m, which is still too shallow to resolve adequately using AEM data.In practice, shallow thaw and sporadic permafrost trends are observed to greater depths in many locations, including inactive or abandoned channels (Jepsen et al., 2013b).To simulate these types of features, the shallow resistivity structure of the 6 m deep hydrostatic lake scenario at 1000 years is manually modified to include three synthetic "channels".These channels are not intended to represent realistic pathways relative to the lake and the hydrologic simulations; they are solely for the purpose of illustrating the ability to resolve shallow resistivity features.Figure 10a shows the three channels in a zoomed-in view of the uppermost portion of the model outside the lake extent.Each channel is 100 m wide but with different depths: 1 m (half the unit 1 thickness), 2 m (full unit 1 thickness), and 3 m (extending into the top of unit 2).Analysis of AEM data simulated for this model, presented as the McMC average model, is shown in Fig. 10b.All three channels are clearly identified, but their thicknesses and resistivity values are overestimated and cannot be distinguished from one another.To explore the possibility of better resolving these shallow features, synthetic electromagnetic data are simulated using the characteristics of a ground-based multi-frequency electromagnetic tool (the GEM-2 instrument) that can be hand carried or towed behind a vehicle and is commonly used for shallow investigations.The McMC average model result for the simulated shallow electromagnetic data is shown in Fig. 10c.An error model with 4 % relative data errors and an absolute error floor of 75 ppm was used for the GEM-2 data.Channel thicknesses and resistivity values are better resolved compared with the AEM result, though the 1 m deep channel near r = 800 m appears both too thick and too resistive.In addition, the shallow electromagnetic data show some sensitivity to the interface at 2 m depth between frozen silty sands and frozen gravels, though the depth of this interface is overestimated due to the limited sensitivity to these very resistive features.

Discussion
Understanding the hydrogeophysical responses to permafrost dynamics under different hydrologic and climatic conditions, and in different geological settings, is important for guiding the interpretation of existing geophysical data sets and also for planning future surveys.Geophysical models are inherently uncertain and ambiguous because of (1) the resolution limitations of any geophysical method and (2) the weak or non-unique relationship between hydrologic properties and geophysical properties.We have presented a general framework for coupling airborne and ground-based electromagnetic predictions to hydrologic simulations of permafrost evolution, including a novel physical property relationship that accounts for the electrical response to changes in lithology, temperature, and ice content, as well as a rigorous analysis of geophysical parameter uncertainty.Although the focus here is on AEM data, other types of electrical or electromagnetic measurements could be readily simulated using the same resistivity model.Future efforts will focus on the simulation of other types of geophysical data (e.g., nuclear magnetic resonance or ground penetrating radar) using the same basic modeling approach.
In the specific examples of lake-talik evolution presented here, which are modeled after the physical setting of the Yukon Flats, Alaska (Minsley et al., 2012), AEM data are shown to be generally capable of resolving large-scale permafrost and geological features (Fig. 5), as well as thermally and hydrologically induced changes in permafrost (Figs. 8  and 9).The Bayesian McMC analysis provides useful details about model resolution and uncertainty that cannot be assessed using traditional inversion methods that produce a single "best" model.A fortuitous aspect of the Yukon Flats model is the fact that the silt layer (unit 3) is relatively conductive compared with the overlying gravels (unit 2), making it a good target for electromagnetic methods.If the order of these layers were reversed, if the base of permafrost were hosted in a relatively resistive lithology, or if the base of permafrost was significantly deeper, AEM data would not likely resolve the overall structure with such good fidelity.In addition, knowledge of the stratigraphy helps to remove the ambiguity between unfrozen gravels and frozen silts, which have similar intermediate resistivity values (Figs. 4 and 5).The methods developed here that use a physical property model to link hydrologic and geophysical properties provide the necessary framework to test other more challenging hydrogeological scenarios.
Two key challenges for the lake-talik scenarios were identified: (1) resolving the details of partial ice saturation beneath the lake during talik formation and (2) resolving nearsurface details associated with shallow thaw.The first challenge is confirmed by Figs. 5 and 8, which show that AEM data cannot resolve the details of partial ice saturation beneath a forming talik.However, there is clearly a change in the overall characteristics of the sub-lake resistivity struc-ture as thaw increases (Fig. 9).One notable feature is the steadily decreasing depth to the top of the low-resistivity unfrozen silt (red) beneath the lake (Fig. 5a-b) as thaw increases, ultimately terminating at the depth of the gravelsilt interface when fully unfrozen conditions exist (Fig. 5c).Measurements of the difference in elevation between the interpreted top of unfrozen silt and the base of nearby frozen gravels were used by Jepsen et al. (2013a) to classify whether or not fully thawed conditions existed beneath lakes in the Yukon Flats AEM survey described by Minsley et al. (2012).The simulations presented here support use of this metric to distinguish full versus partial thaw beneath lakes.However, without the presence of a lithological boundary, the shallowing base of permafrost associated with talik development beneath lakes would be much more difficult to distinguish.Finally, it is important to note that resistivity is sensitive primarily to unfrozen water content and that significant unfrozen water can remain in relatively warm permafrost that is near 0 • C, particularly in fine-grained sediments.Resistivityderived estimates of talik boundaries defined by water content may therefore differ from the thermal boundary defined at 0 • C.
The second challenge, to resolve near-surface details associated with supra-permafrost thaw, is addressed in Fig. 10.For the scenarios considered here, AEM data can identify shallow thaw features but have difficulty in discriminating their specific details.There are many combinations of resistivity and thickness that produce the same electromagnetic response; therefore, without additional information it is not possible to uniquely characterize both thaw depth and resistivity.Ground-based electromagnetic data show improved sensitivity to the shallow channels and also limited sensitivity to the interface between resistive frozen gravels and frozen silty sands (Fig. 5).By restricting the possible values of resistivity and/or thickness for one or more layers based on prior assumptions, Dafflon et al. (2013) showed that improved estimates of active layer and permafrost properties can be obtained.The quality of these estimates, of course, depends on the accuracy of prior constraints used.In many instances, it may be possible to auger into this shallow layer to provide direct observations that can be used as constraints.This approach could be readily applied to the ensemble of McMC models.For example, if the resistivity of the channels in Fig. 10a were known, the thickness of the channels could be estimated more accurately by selecting only the set of McMC models with channel resistivity close to the true value, thereby removing some of the ambiguity due to equivalences between layer resistivity and thickness.
AEM data are most likely to be useful for baseline characterization of subsurface properties as opposed to monitoring changes in permafrost.Although there are some cases of rapid change associated with near-surface freeze-thaw processes (Koch et al., 2013), and the case of catastrophic loss of ice in the gaining lake scenario (Fig. 9b), that may be of interest, large-scale changes in permafrost generally occur over much longer time periods than is practical for repeat AEM surveys.One exception could be related to infrastructure projects such as water reservoirs or mine tailing impoundments behind dams, where AEM could be useful for baseline characterization and repeat monitoring of the impact caused by human-induced permafrost change.Geophysical modeling, thermophysical hydrologic modeling, and field observations create a synergy that provides greater insight than any individual approach and can be useful for future characterization of coupled permafrost and hydrologic processes.

Summary
Analysis of AEM surveys provide a means for remotely detecting subsurface electrical resistivity associated with the co-evolution of permafrost and hydrologic systems over areas relevant to catchment-scale and larger processes.Coupled hydrogeophysical simulations using a novel physical property relationship that accounts for the effects of lithology, ice saturation, and temperature on electrical resistivity provide a systematic framework for exploring the geophysical response to various scenarios of permafrost evolution under different hydrological forcing.This modeling approach provides a means of robustly testing the interpretation of AEM data given the paucity of deep boreholes and other ground truth data that are needed to characterize subsurface permafrost.A robust uncertainty analysis of the geophysical simulations provides important new quantitative information about the types of features that can be resolved using AEM data given the inherent resolution limitations of geophysical measurements and ambiguities in the physical property relationships.In the scenarios considered here, we have shown that large-scale geologic and permafrost structure is accurately estimated.Sublacustrine thaw can also be identified, but the specific geometry of partial ice saturation beneath lakes can be poorly resolved by AEM data.Understanding the geophysical response to known simulations is helpful both for guiding the interpretation of existing AEM data and to plan future surveys and other ground-based data acquisition efforts.

Figure 1 .
Figure 1.Axis-symmetric model geometry indicating different lithologic units and simulated lake depths/extents.
Figure2.Bulk resistivity as a function of ice saturation using the physical properties defined for each of the lithologic units described in Table1.

Figure 3 .
Figure 3. SUTRA model outputs and geophysical transformations from the 6 m gaining lake simulation at 240 years.Ice saturation (a) and temperature (b) are converted to predictions of bulk resistivity (c).Variability in resistivity as a function of temperature is indicated in (d) for lithologic units 1-3.

Figure 4 .
Figure 4. Synthetic bulk resistivity images under hydrostatic flow and current climate conditions.Lake depths of 3 m (a-c), 6 m (d-f), 9 m (g-i), and 12 m (j-l) are illustrated at simulation times 100, 240, and 1000 years.

Figure 5 .
Figure 5. Mean resistivity model extracted from McMC ensembles.Results are shown for the 6 m deep hydrostatic lake scenario outputs at (a) 100 years, (b) 240 years, and (c) 1000 years.

Figure 6 .
Figure 6.Performance of geophysical parameter estimation in recovering true parameter values.(a) True versus McMC-estimated resistivity values for the hydrostatic 6 m deep lake scenario at simulation time 1000 years compared with the frequency distribution of true resistivity values (b).Estimated resistivity values generally fall along the dashed 1 : 1 line in (a), with the exception of underprediction of the resistive frozen gravels, overprediction of the thin surficial frozen sand, and some overprediction of the frozen silt where it is in contact with frozen gravel.

Figure 7 .
Figure 7. McMC-estimated resistivity posterior distributions within frozen and unfrozen unit 2 gravels (a) and frozen and unfrozen unit 3 silts (b) for the hydrostatic 6 m deep lake scenario at years.Unfrozen resistivity distributions are extracted beneath the center of the lake (r = 0) at depths of 15 and 50 m for the gravels and silts, respectively.Frozen distributions are extracted at the same depths but at r = 750 m.The upper x-axis labels indicate approximate ice saturation based on the lithology-dependent ice saturation versus resistivity curves shown in Fig. 2.

Figure 8 .
Figure 8. Resistivity probability distributions for the hydrostatic 6 m deep lake scenario at simulation times 100 years (a-b), 240 years (c-d), and 1000 years (e-f).Shading in each image represents the probability distribution at depths of 15 m (a, c, e) and 50 m (b, d, f) from the lake center (r = 0 m) to the edge of the model (r = 1800 m).Dashed lines indicate the true resistivity values.Ice saturation is displayed on the right axis of each image and is defined empirically for each lithology using the relationships in Fig. 2.

Figure 9 .
Figure 9. Change in ice saturation and resistivity as a function of time.Results are shown for the 6 m deep lake hydrostatic and gaining lake scenarios within (a) the gravel layer, unit #2, at a depth of 15 m and (b) the silt layer, unit 3, at a depth of 50 m.

Figure 10 .
Figure 10.Comparison of airborne and ground-based measurements for recovering shallow thaw features.(a) True shallow resistivity structure extracted from the hydrostatic 6 m deep lake scenario at a simulation time of 1000 years, shown outside of the lake extent (distance > 500 m).Three shallow low-resistivity channels with thicknesses of 1, 2, and 3 m were added to the resistivity model to provide added contrast.McMC-derived results using simulated AEM data (b) and ground-based electromagnetic data (c) illustrate the capability of these systems to image shallow features.

Table 1 .
Description of geologic units and physical properties used in numerical simulations.Entries separated by commas represent parameters with different values for each of the lithologic units.