Interactive comment on “ Sensitivity of airborne geophysical data to sublacustrine permafrost thaw ”

Introduction Conclusions References


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;Figures Back Close Full  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 characteriz-Introduction

Conclusions References
Tables Figures

Back Close
Full ing 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., 2012a), 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 (sNMR) 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 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 coupled thermal-hydrologic simulations of Wellman et al. (2013) predict the evolution of lake taliks (unfrozen sublacustrine areas in permafrost regions) 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) simulation of synthetic geophysical data that would be expected for various permafrost hydrologic conditions; 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., 2012a); a lake that is also the basis for the lake simulations discussed by Wellman et al. (2013).Introduction

Conclusions References
Tables Figures

Back Close
Full  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 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 m 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., 2012a;Williams, 1962), with defined hydrologic and geophysical parameters for each layer summarized in Table 1.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 tex-Introduction

Conclusions References
Tables Figures

Back Close
Full ture, 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 i th species, respectively.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 conduc-Introduction

Conclusions References
Tables Figures

Back Close
Full tivity) 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, 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, 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

Conclusions References
Tables Figures

Back Close
Full 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 airborne electromagnetic (AEM) data are simulated for each snapshot of predicted bulk resistivity values using nominal system parameters based on the Fugro RESOLVE1 frequency-domain AEM system that was used in the Yukon Flats survey (Minsley et al., 2012a).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 (VCX) coil pair with 9 m separation that operates at 3.260 kHz.Oscillating currents and associated magnetic fields in 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 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 equations described in Minsley (2011).Introduction

Conclusions References
Tables Figures

Back Close
Full 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) 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 5-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 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 0 and 1000 years output 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 hand-carried 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 EM data (Minsley, 2011) to explore the ability of simulated AEM data to recover the true distribution of subsurface resistivity values throughout the 1000 year lake talik simulations.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 sampled 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

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 SUTRA 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 vs. 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 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 zero degrees, with a rapid increase in resistivity for temperatures below zero degrees.This result is very similar to the temperature-resistivity relationships illustrated by Hoekstra (1975, Fig. 1), lending confidence to our physical property definitions described earlier.Above zero degrees, the slight decrease in resistivity is due to the temperaturedependence of fluid resistivity.The rapid increase in resistivity below zero degrees is primarily caused by reductions in effective porosity due to increasing ice saturation, though changes in surface conductivity and salinity at increasing ice saturation are also contributing factors.Below −1 • C, the change in resistivity values as a function of temperature rapidly decreases.This is an artifact caused by the imposed temperatureice saturation relationship defined in SUTRA that, for these examples, enforces 99 % ice saturation at −1 • C. It is more likely that ice saturation continues to increase asymptotically over a larger range of temperatures below zero degrees, with corresponding increases in electrical resistivity.However, because AEM methods are limited in 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
Geophysical data (not shown) are simulated for each of the electrical resistivity models (e.g.Fig. 4) using the methods described in Sect. the average resistivity model as a function of depth is calculated from the McMC ensemble of 100 000 plausible models.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) vs. 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 vs. 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 EM 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.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, whose resistivity 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 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 (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 mdeep 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 causes 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.In addition, the shallow EM 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 over-estimated 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 datasets 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 geophysical 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 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., 2012b), AEM data are shown to be generally capable of resolving large-scale permafrost and geological

Conclusions References
Tables Figures

Back Close
Full features (Fig. 5), as well as thermally and hydrologically induced changes in permafrost over time (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, or if the base of permafrost were hosted in a relatively resistive lithology, 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 near-surface 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 structure 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 and b) as thaw increases, ultimately terminating at the depth of the gravel-silt 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. (2012a).The simulations presented here support use of this metric to distinguish full vs. 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.

Conclusions References
Tables Figures

Back Close
Full The second challenge, to resolve near-surface details associated with suprapermafrost 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 EM response; therefore, without additional information it is not possible to uniquely characterize both thaw depth and resistivity.Ground-based EM 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) or in the case of catastrophic loss of ice in the gaining lake scenario (Fig. 9b) that may be of interest, large-scale changes in permafrost occur over much longer time periods than is practical for repeat AEM surveys.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.

Conclusions References
Tables Figures

Back Close
Full

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 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 quantitative information about what types of features 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 also to plan future surveys and other ground-based data acquisition efforts.defined for each of the lithologic units described in Table 1.746 747 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3 Results Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2.3.The simulated data are then used to recover estimates of the original resistivity values according to the approach outlined in Sect.2.4.Resistivity parameter estimation results for the 6 m-deep hydrostatic lake scenario (Fig. 4d-f) are shown in Fig. 5.At each location along the profile, Discussion Paper | Discussion Paper | Discussion Paper | 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 ice-saturation values, translated from the ice vs. 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: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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, are 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 Discussion Paper | Discussion Paper | Discussion Paper |features, synthetic EM data are simulated using the characteristics of a ground-based multi-frequency EM 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 EM data is shown in Fig.10c.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.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 6 .Figure 7 .
Figure 6.Performance of geophysical parameter estimation in recovering true parameter values.(a) True vs. 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 exceptions being under-prediction of the resistive frozen gravels, over-prediction of the thin surficial frozen sand, and some over-prediction of the frozen silt where it is in contact with frozen gravel.

Figure 7 .Figure 8 .
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 1000 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 axes labels indicate approximate ice saturation based on the lithology-dependent ice saturation vs. resistivity curves shown in Fig. 2.

Figure 8 .Figure 9 .Figure 9 .Figure 10 .
Figure 8. Resistivity probability distributions for the hydrostatic 6 m-deep lake scenario at simulation times 100 years (a and b), 240 years (c and d), and 1000 years (e and 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 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 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 EM data (c) illustrate the capability of these systems to image shallow features.