Mapping Ice Buried by the 1875 and 1961 Tephra of Askja Volcano, Northern Iceland Using Ground‐Penetrating Radar: Implications for Askja Caldera as a Geophysical Testbed for In Situ Resource Utilization

Eruptions of the Askja Volcano in Northern Iceland in 1875 and 1961 blanketed the caldera with rhyolitic and basaltic tephra deposits, respectively, which preserved layers of seasonal snowpack as massive ice. Askja serves as an operational and geophysical analog to test ground‐penetrating radar field and analysis techniques for in situ resource utilization objectives relevant to the martian and lunar environments. We conducted ground‐penetrating radar surveys at center frequencies of 200, 400, and 900 MHz to map the thickness and extent of tephra deposits and underlying massive ice at three caldera sites. We identified up to 1 m of tephra preserving up to 4.4 m of massive ice. We measured the real dielectric permittivity of the overlying tephra and the total attenuation at each frequency of the tephra and ice. A key objective of our investigation was to determine if attenuation (or loss) could be used as an additional diagnostic signature of massive ice preserved at depth when compared to ice‐free stratigraphy. Loss rates of the ice‐rich subsurface decrease with increasing ice thickness relative to the overburden, which may constitute a possible signature. Attenuation also increased with increasing frequency. The tephra, ice, and other volcanic deposits at each of our three caldera sites and the ice‐free, pumice‐mantled 1961 Vikrahraun lava flow exhibited consistently low loss rates at all frequencies. This result highlights the ambiguity associated with identifying the unique signature of ice within low‐loss stratigraphies, a possible challenge for its identification in the martian or lunar subsurface using radar.


Introduction
Ground-penetrating radar (GPR) surveying is widely employed in periglacial and glacial settings on Earth to characterize the subsurface physical properties of ice, soil, and sediment in the upper tens of meters of the subsurface.GPR has been identified as a non-invasive and easily deployable tool to investigate the subsurface and characteristics of regolith and ice-containing sites on the Moon and Mars (e.g., Boisson et al., 2011;Grimm et al., 2006;Hamran et al., 2020Hamran et al., , 2022;;Heggy, Clifford, Grimm, Dinwiddie, Stamatakos, et al., 2006;Heggy, Clifford, Grimm, Dinwiddie, Wyrick, et al., 2006;Lai et al., 2019;Li et al., 2020Li et al., , 2022;;LWIMS Report, 2020;MEPAG ICE-SAG Report, 2019;Richardson et al., 2020;Shoemaker et al., 2022 and references therein).Orbital radar sounders were the first systems used to identify and confirm widespread, buried ice deposits across the midlatitudes of Mars (e.g., Bramson et al., 2015;Dundas et al., 2018;Holt et al., 2008;Morgan et al., 2021;Picardi et al., 2004;Seu et al., 2007;Stuurman et al., 2016).GPR systems have now been successfully deployed on rovers on the surfaces of Mars and the Moon (Fang et al., 2014;Hamran et al., 2020Hamran et al., , 2022;;Li et al., 2020Li et al., , 2022)).There is significant interest in including a GPR system as part of a future in situ resource utilization (ISRU) campaign to a polar cold trap on the Moon, similar in scope to the upcoming VIPER mission (Volatiles Investigating Polar Exploration Rover) (Colaprete, 2021;LWIMS Report, 2020;Shoemaker et al., 2022).
Terrestrial analog field sites are critical to test GPR methods to successfully characterize the properties of subsurface regolith and ice for future ISRU campaigns (Boisson et al., 2011;Brandt et al., 2007;Campbell et al., 2018;Dinwiddie et al., 2005;Grimm et al., 2006;Heggy, Clifford, Grimm, Dinwiddie, Stamatakos, et al., 2006, Heggy, Clifford, Grimm, Dinwiddie, Wyrick, et al., 2006).At an elevation of >1 km, Askja is in a region of discontinuous permafrost (Czekirda et al., 2019;Etzelmüller et al., 2007Etzelmüller et al., , 2020) ) that has been uniquely influenced by its regional volcanic activity and historical deposits of tephra (Kellerer-Pirklbauer et al., 2007).Askja provides a geophysical testbed to probe shallow (0-10 m) ice deposits similar in thickness and burial to those on Mars or potentially the Moon.Analytical and field deployment methods can be readily developed there to achieve future science and ISRU objectives for terrestrial bodies (Cannon & Britt, 2020;Ellery, 2020;Starr & Muscatello, 2020).A current challenge in utilizing GPR for coordinated resource campaigns is the ambiguous results it can yield in subsurface stratigraphy.This has been demonstrated for layers such as water ice and a lowdensity regolith, which share similar dielectric properties (Boisson et al., 2011).Different depositional processes may also result in similar morphologies in radar returns ("radar facies"), which can complicate interpretations that guide real-time science objectives of rover surface missions (e.g., Hamran et al., 2022).Our investigation seeks to determine if the total losses to the radar signal can be used as an additional diagnostic signature of the presence of subsurface water ice when compared to ice-free locations, given several meters of buried ice at Askja.To accomplish this, we conducted the first multi-frequency GPR mapping campaign to characterize the thickness, extent, and dielectric properties of the massive ice and pyroclastic deposits at the Askja Volcano located in the Northern Icelandic Highlands (Figure 1a).

Askja Eruptions and Caldera Site Description
The Askja central volcano is located in the Northern Volcanic Zone of Iceland (Figure 1a,inset).Askja was the source of a series of recent explosive and effusive eruptions in 1875, 1921-1922, 1929, 1931, and 1961 in addition to many other earlier Holocene eruptions (Annertz et al., 1985;Carey et al., 2008aCarey et al., , 2008b;;Graettinger et al., 2013;Hartley et al., 2016).The youngest and current caldera is Öskjuvatn (Figure 1a), now occupied by a lake, formed by an explosive phreatoplinian eruption in March 1875 (Carey et al., 2009(Carey et al., , 2010;;Graettinger et al., 2013;Sparks et al., 1981).The March 1875 eruption produced 0.33 km 3 dense rock equivalent of tephra throughout eastern Iceland and parts of Scandinavia and Germany (Carey et al., 2009).The caldera tephra deposit is a coarse, sub-angular, rhyolitic pumice (referred to as Unit D; Sparks et al., 1981) along with an underlying unit of gray/brown massive ash fall and/or intracaldera pyroclastic surge deposits (broadly referred to as Unit C; Carey et al., 2010;Sparks et al., 1981).In October-November 1961, vents near the northern rim of the caldera deposited a black-brown basaltic lapilli and some ash throughout the caldera and the Vikrahraun lava flow was emplaced eastward from the rim (Blasizzo et al., 2022;Thorarinsson & Sigvaldason, 1962).These two tephra deposits each blanketed and preserved a layer of seasonal snowpack that was later densified into massive ice (an extensive layer comprised mostly of lithic-poor ice) (Carey et al., 2009;Helgason, 2000).This permafrost aggradation process initiated via snowpack burial and preservation by volcanic tephra has been observed elsewhere in Iceland in its early stages at the Hekla volcano after its 2000 eruption and possibly in its much later stages at Öraefajökull preserved by ash from its eruption in 1362 (Helgason, 2000;Kellerer-Pirklbauer et al., 2007).Various thermokarst and permafrost landforms are observed across the Askja caldera and flanks that indicate the presence of ice at depth.We focus on three sites within the caldera that are blanketed by these tephra deposits (Figure 1b, white boxes).We chose these sites because they possess morphologic evidence that buried ice is present and represent the major tephra deposits and near-surface stratigraphy observed within this region of the caldera.Figures 1c-1e show each of these three sites and associated GPR surveys taken in 2019 and 2021 (yellow lines).From north to south, Site 1 (Figures 1c and 2a) possesses massive ice buried by the 1961 tephra deposit, Site 2 (Figures 1d and 2b) captures the 1961 tephra beginning to transition to the 1875 tephra, and Site 3 (Figures 1e and 2c) is primarily blanketed by the 1875 tephra with 1961 tephra mostly confined to topographic depressions.These three sites exhibit a range of tephra overburden types and clast sizes, ice thickness, and permafrost landforms, making this region of Askja an excellent site to test GPR field and data analysis techniques.

Field Methods
We conducted two GPR field campaigns within the Askja caldera in July-August 2019 and in August 2021 on rain-free days.At each site during both campaigns, relatively dry conditions and porous tephra kept moisture levels in the subsurface to a minimum, except for a limited zone of elevated moisture observed at the base of the tephra layer where it contacts the ice table.We used a Geophysical Survey Systems Inc. (GSSI) SIR 4000 GPR system with three shielded antennas operating at center frequencies of 200, 400, and 900 MHz for surveys.Crosssectional images of the subsurface (radargrams) were collected as a series of individual traces along the length of a traverse across the surface.Traces were collected in distance mode by an attached odometer wheel where the profile of the subsurface was sampled at 100 traces/meter.Radargrams are displayed as two-way travel times (or depth) of returned radar wave amplitudes as a function of along-traverse distance (Figure 3).Strong contrasts in dielectric constant or the real relative permittivity between subsurface materials will appear as a distinct bright boundary or reflector at depth.These reflections allowed us to map the horizontal and vertical extent of tephra and ice layers across the Askja caldera.
Position and altitude of the GPR unit were controlled using a combination of differential GPS (dGPS) and contextual high-resolution aerial surveys of the caldera floor using a small uncrewed aerial system (UAS) DJI Mavic Pro.Traverses were georegistered and terrain-corrected using a Trimble Geo7x Handheld GPS attached to the GPR unit with a tens-of-centimeter-to-meter position accuracy.UAS surveys took place in 2019 and 2021 to provide detailed knowledge of surface lithology, water and ice, terrain height, and additional positional control of GPR transects.UAS images were used to produce orthomosaics (Figures 1c-1e) and We mapped massive ice buried by tephra using ground-penetrating radar (GPR) at center frequencies of 200, 400, and 900 MHz.Panels a, b, and c are ice thickness maps summarized for the 400 MHz observations.Massive ice was observed to be thickest at Site 2, panel (b), reaching depths of 4.4 m.Ice thicknesses <10 cm are reported as zero on these maps.Radargrams collected at each center frequency at each of the three sites are summarized in Figure 3. Trench and borehole locations are marked in red.(a) Site 1 massive ice deposits are thinner in the west and thicker toward the hiking trail in the east, reaching a maximum depth of 2.8 m.The tephra transitions from more loosely packed ripples or aeolian bedforms to a smoother, compact section after point A. Internal layering was observed within the less compact tephra.(b) Site 2 has more ablation features than Site 1 (e.g., tension cracks and depressions).Massive ice is likely resting on an older ash or tephra deposited prior to 1961.Massive ice deposits are also thickest at this site, reaching a maximum of 4.4 m moving south along the hiking trail.(c) Massive ice at Site 3 tends to be thicker, on average, beneath the 1875 pumice.Obvious thermokarsts are scattered throughout this region.A meltwater channel generated by seasonal melting of snow cuts through GPR traverses taken on either side.Trenches taken on the east and west banks of this channel revealed ice buried by tephra from the 1875 and 1961 eruptions.Images of the interplay of these deposits and massive ice along with labeled trenches and boreholes from this figure are shown in Figure 4. corresponding digital elevation models of our survey sites.The 2019 mosaic and elevation model data products were produced at a spatial resolution <3 cm/pixel in AgiSoft Metashape software; 2021 data products were produced at <6 cm/pixel.UAS data products were georegistered to a horizontal and vertical precision <5 cm/pixel with a post-processing kinematic dGPS survey conducted at stationary markers laid out in the UAS Survey area.Additionally, we drilled boreholes and dug trenches along collected GPR traverses to collect observations of tephra stratigraphy variations across the caldera and to verify subsurface reflectors identified in radargrams, locate the top of the ice table at survey sites, and note changes in moisture conditions with depth.

Quantifying Permittivity and EM Wave Velocity
The permittivity of geologic materials is an intrinsic property expressed as a complex number: ε* = ε′ iε″.The ratio of the imaginary component, ε″, describing signal losses, and the real component, ε′, is the loss tangent, tanδ.This quantity is, in turn, related to the attenuation of the radar signal (see Section 3.4).Relative permittivity is given by ε r = ε/ε 0 where ε 0 is the permittivity of free space.The real component of the relative permittivity, ε′ r , (or dielectric constant; hereafter "permittivity" for brevity) is related to the radar wave propagation through a medium, which is given by where v is the bulk velocity of the radar wave through the medium and c is the speed of light in vacuum.ε′ r of subsurface layers was estimated at different sites using two methods.First, we estimated the velocity and calculated permittivity (Equation 1) from trench measurements of tephra thickness and the one-way travel times from radar picks at the tephra-ice interface.Second, we estimated bulk velocity by fitting theoretical hyperbolas to hyperbolic forms that are generated by the motion of the GPR system toward, over, and away from an embedded subsurface scatterer.We use the GSSI RADAN 7 processing and analysis software to fit hyperbolas where possible and report estimates of bulk velocity and permittivity for those sections of the subsurface above the hyperbola sources.

Mapping Ice and Tephra Thicknesses
We produced maps of ice and tephra thickness at each site by picking radar reflectors within collected radargrams.We applied a time-zero correction and an exponential gain to the radargrams prior to picking.This gain was applied only for visualization purposes and removed for any calculations.We used the open source Radar Analysis Graphical Utility (RAGU) picking software (Tober & Christoffersen, 2020) to generate our reflector inventory.From each picked reflector, we exported relevant quantities to an ESRI shapefile, including radar amplitudes, latitude and longitude, two-way travel times, elevation, and layer thickness.Layer thicknesses were calculated using permittivity values estimated using the two methods described in Section 3.2, with a trench or borehole measurements preferred over hyperbola fits.Each shapefile was then input into GIS software where it was overlaid onto the caldera orthoimages.Examples of these ice thickness maps are shown in Figure 2 and Figure S1 in Supporting Information S1.The top of the ice table is fairly constant, measured from trenching to be at 52 cm depth beneath the 1961 tephra.From the ice table position, an average permittivity of 2.91 is applied to place the ice table at its observed depth in the images.(b) Site 2 radargrams correspond to panel (b) of Figure 2. Section B to B' is a 900 MHz observation with relatively thin ice (∼50 cm) preserved by ∼33 cm of 1961 tephra taken in 2021.This radargram was depth-corrected using a derived tephra permittivity of 7. Section C to C′ is repeat 400 and 200 MHz observations taken in 2019 where up to 65 cm of massive ice is preserved by ∼50 cm of 1961 tephra.These radargrams were depth-corrected using a permittivity of 2.59.(c) Site 3 radargrams correspond to panel (c) of Figure 2. Radargrams from D to D′ were taken upslope.The cm-to-dm-sized clasts of pumice from 1875 are visibly more scattering than the sub-wavelength 1961 tephra.Massive ice reaches a maximum thickness for Site 3 along this section at 2.8 m.The ice table is estimated to be at ∼62 cm from a borehole 29 m along-traverse.These radargrams were depth-corrected using an average tephra permittivity of 4.6.All radargrams were processed in the GSSI RADAN 7 software using a vertical infinite impulse response filter, a five-point exponential gain, and corrected for variations in surface topography.In panel (c), ∼15 cm black/brown lapilli and coarse ash from 1961 give way to 15 cm of moist brown ash with mm-to-cm-sized blocks of buff colored 1875 pumice.The ice table rests at 30 cm depth.T2 in panel (d) was dug on a slope of the eastern meltwater channel bank in Figure 2c, making observed depths larger than vertical.∼25 cm coarse, mm-to-cm-size clasts of black/brown 1961 tephra and gray/brown ash transitions to ∼15 cm of pore ice and 70 cm of massive ice.∼53 cm of an ice-cemented, pumice-supported coarse-to-fine ash matrix from 1875 extends to the floor of the channel, which is comprised of a 50/50 mixture of pumice clasts and ice from 1875.(e) Borehole B7 (Figure 2c) was taken at 29 m along radargrams in Figure 3c.This retrieval (arranged in discontinuous but depth-ordered sections) shows the transition from ice-cemented 1875 pumice and ash to massive ice from bottom to top in the image.The ice table was measured to be at 62 cm depth.The trench was taken along the eastern meltwater channel bank near these GPR traverses.∼40 cm of cm-sized 1875 pumice with ash is burying ∼3 m of massive ice at Site 3. (f) This trench was dug along the western bank of the meltwater channel at Site 3 where 36 cm of 1961 basaltic lapilli was covering ∼60 cm of massive ice that was resting on moist, brown 1875 ash deposits with some mm-to-cm-sized clasts of buff colored 1875 pumice.(g) GPR traverses across the toe of the 1875 pumice-mantled 1961 Vikrahraun lava flow found ∼1 m of pumice rests atop the basalt flow.1875 tephra was finer-grained at this site (mm-sized grains and fine-grained ash).Stratigraphy is ice-free.

Quantifying Losses to the Radar Signal
Total losses to the radar signal are sourced from intrinsic attenuation (i.e., absorption), scattering, and the geometric spreading of the transmitted wavefront as it travels through the subsurface and back to the receiver.Losses are directly proportional to distance, R, and exhibit a semilogarithmic decay in amplitude with depth (or travel time) given by e 2αR , for uniform layer thicknesses and constant reflection coefficients, where α is the attenuation coefficient (Grimm et al., 2006).This behavior can also arise from an isotropic distribution of scatterers (Grimm et al., 2006).Similar behavior has been observed in seismology (e.g., Farrokhi et al., 2016;Jin et al., 1994) and in both cases, this constant decay is referred to as Q where Q 1 = tan δ = αλ/π, where tan δ is loss tangent and λ is the wavelength.To quantify the losses to the radar signal, we employ methods similar to Grimm et al. (2006) and Boisson et al. (2011) and fit sections of averaged amplitudes that exhibit a semilogarithmic decay with depth.To estimate the total loss and Q from the amplitudes at each GPR antenna frequency, we first applied a series of processing steps to the GPR data, including (in order): (a) time-zero correction, (b) horizontal background filter to remove sources of coherent backscatter, (c) Hilbert transform to obtain the magnitude of collected traces and to further reduce signal variation, and (d) amplitude-normalization using the peak of the direct wave.From the processed data, we then calculated the average of traces over a segment of the radargram where layers were constant in time-delay and maintained uniform thickness.
Total attenuation is estimated from trace averages corrected for ground losses such as the geometric spreading of the wavefront and the backscatter cross sections of the reflecting targets.These losses are described by the radar equation (Annan & Davis, 1977;Skolnik, 2008), the ratio of the reflected power to the transmitted power: where G is antenna gain, λ is the radar wavelength, ξ is the backscatter cross-section, α is the intrinsic attenuation coefficient, and R is the distance to the reflecting target.After correcting the relative amplitudes for these effects, total attenuation can be isolated from the average trace (Annan & Davis, 1977;Boisson et al., 2011;Grimm et al., 2006;Scabbia & Heggy, 2018).
Three models are considered for the backscatter cross-section, ξ: (a) a smooth, planar reflector where ξ = πR 2 Γ where Γ is the power reflection coefficient, therefore, P R /P T ∝ 1/R 2 , (b) GPR returns are integrated over the diameter of the first Fresnel zone (
We applied all three models to correct the trace averages; however, the rough reflector and Rayleigh scatterer models overcorrected the data, which generated an unrealistic positive slope.We therefore applied a modeldependent gain function of the form 1/R 2 to the average trace, where R is the depth in meters.Round-trip travel times are converted to depth using a three-layer permittivity model.Permittivity of the overlying tephra is fairly well constrained from trenches and at the base of that layer (see Sections 3.2 and 4.2).We assume a permittivity of 3.2 for ice at temperatures at or around 0°C for microwave frequencies (Johari, 1976;Matsuoka et al., 1997).Our assumption of the permittivity of the deepest layer then depends on our observations at each site, which is further discussed in Section 4.3 (see Table S3 for layer permittivity).One-way loss rates in dB/m are estimated from linear least squares fits to the portion of the corrected average trace exhibiting an exponential decay with depth, from the position of the first positive peak of the corrected average trace to the noise floor (see Figures 5a-5d).We used the slope from our fits to estimate loss tangent where tan δ 1 = Q.We used this total attenuation to test whether the presence of massive ice at depth would result in a lower loss when compared to icefree stratigraphy.
Estimates of loss rates, Q, and loss tangent were made using open-source Python code developed by Shoemaker (2023).This code depends on the open-source Python software readgssi (Nesbitt et al., 2022) and functions from the open-source Python software RAGU (Tober & Christoffersen, 2020) to read in and auto-detect the peak of the direct wave from each collected GSSI radargram prior to averaging traces, correcting the total loss curve for spreading and scattering effects, and fitting.There is a significant spread in the losses associated with each deposit at each frequency.We report the slopes of successful fits along with other parameters in Table S3.

Distribution of Ice and Tephra at Surveyed Sites
We examined 57 radargrams collected across the three caldera sites at central frequencies of 200, 400, and 900 MHz (Figure 1).Dominant tephra layers at Sites 1-3 are the 1961 tephra and 1875 tephra Units C and D (Carey et al., 2010;Sparks et al., 1981).Unit D is a Plinian fall deposit of well-sorted lapilli.We observed some larger cm-to-dm-sized pumice clasts at Site 3. Unit D thicknesses of up to 200 cm proximal to Öskjuvatn are estimated by Carey et al. (2010).Unit C underlies D and is a gray/brown massive ash fall up to ∼20-25 cm thick in Site 3 and potentially present at Sites 1 and 2 at depth.Intracaldera pyroclastic surge deposits comprise subunits of C but are confined to regions proximal to Öskjuvatn and not present at our sites of interest.
Reflections associated with the top of subsurface massive ice (and base of the tephra layer(s)) were identified in GPR data at each of the sites.We summarize the maximum tephra and massive ice thicknesses and their averages for all sites at 400 MHz in Table 1.
Corresponding ice thickness estimates for each of the three sites at 400 MHz (for clarity) are shown in Figure 2. A summary of ice thickness for the collected radargrams at each site and at all frequencies can be found in Figure S1 in Supporting Information S1.Ice thickness could not be estimated for all GPR traverses.In limited cases, it was difficult to identify the ice layer at depth due to scattering in the capping tephra layer or poor coupling to the surface by the antenna.We excluded 900 MHz data at Site 3 from ice thickness estimates because of this.In other cases, the lower resolution at 200 MHz prevented the identification of thin layers of tephra.A discussion of observations at individual sites follows this section.

Site 1 Observations
Site 1 contains GPR observations primarily taken during the 2021 field season.The 1961 basaltic tephra is visible at the surface as black/brown lapilli with some gray/brown ash mixed in.A series of long traverses (up to 206 m) were repeated at 200, 400, and 900 MHz over the 1961 basaltic tephra deposit (Figures 1c, 2a, and 4a).From visual observations taken on the surface, this tephra deposit is uniform in physical properties such as appearance, grain size, porosity, and compaction along each survey line.Trenches were dug along the GPR traverse lines, with locations shown in Figure 2. At trench T10 (Figure 4a), the black/brown basaltic 1961 lapilli extends to the top of the ice table at ∼40 cm depth.The subsurface was particularly dry at this site.An ∼5-10 cm thick zone of moisture was observed at the base of the tephra layer in the trench.Isopach maps from Carey et al. (2010) indicate that 1875 Unit D and potentially C may be present at depth.We did not observe any tephras in the trenches at Site 1 that matched the descriptions of these units.It is therefore likely that, if present, they were emplaced prior to snowpack that was later preserved by the 1961 tephra and are stratigraphically lower at this site.
Tephra in the southernmost portion of Site 1 is slightly finer-grained, having been remobilized to form ripples or aeolian bedforms.In radargrams, this area of ripples on the surface corresponds to a series of closely spaced shallow reflectors within the 1961 tephra deposit interpreted to result from internal layering and moisture differences.These reflectors are confined to regions with surface ripples at the start and end of the longest traverse in Figure 2a.The ice deposit is thickest (∼1.4-2.8 m thick) in the south (beneath these bedforms) and east (toward the hiking trail).Indeed, the reflectors identified at the top and base of the massive ice deposit are seen to extend beneath these regions of surface ripples.The transition from these bedforms to a smoother, more compact tephra is visible in Figure 2a just before point A. Massive ice layers at Site 1 are discontinuous in radargrams (Figure 3a), and are thickest in troughs within what is interpreted as an older underlying lava flow.Ice pinches out toward structures interpreted to be buried lava festoons (arcuate surface ridges at m-to-tens-of-m spacing formed by increases in viscosity as lava cools while the interior continues to flow).Trenching revealed the ice table at fairly consistent depths between 40 and 52 cm, except where it was disrupted by lava festoons.Some lava festoons are exposed along the survey line (e.g., near point A) as rocky outcrops at the surface, where they also generated clutter in the GPR data.It may be that 1875 tephras (Units C and D) are present beneath the ice table and overlying these buried lava flows.We do not see multiple reflections associated with these units if they are present; they may be thin enough to be at the resolution limit of the radar antennas used.The mm-to-cm-sized clasts of lapilli generated little scattering visible in the radargrams, as compared to Site 3 (see Section 4.1.3)where most of the ice is buried by the cm-to-dm-sized 1875 pumice.

Site 2 Observations
1961 basaltic tephra covers the surface of Site 2 with some windblown 1875 pumice fragments scattered across the surface (Figures 1d and 2b).GPR surveys were taken across regions almost exclusively covered by 1961 tephra.Ablation and melt-related features are more widespread here than at Site 1, observed as collapsed depressions, tension cracks, and hummocks along with pooled surface water.A borehole taken at this site east of the hiking trail revealed this meltwater was present in the subsurface as well (red point, Figure 2b).Many of these depressions appear to be concentrated or initiated near the boundaries between subsurface ice and abutting lava outcrops.Surveys were taken along the trail and in a region west of the trail where melting has generated several "islands" of massive ice at depth that are surrounded by collapse features and small hummocks (Figure 2b, B to B′).
Massive ice layers are fairly continuous in the radargrams (Figure 3b) but are disrupted by ablation and where collapse features are evident at the surface.The tephra cover is consistently thinner at this site, which could be the reason for the increase in ablation features compared to Site 1. Ice deposits are thicker in the east toward the hiking trail, approaching 5 m and are as thin as 9 cm in the western portions.The ice table was at depths of 33-50 cm in trenches, and is well resolved at each frequency.An ∼10-12 cm thick zone of elevated moisture at the base of the overlying tephra layer was observed in the trenches here (Figure 4b).The massive ice layer in each radargram shows very little scattering or internal reflectors, implying low lithic content and/or dispersed, entrained sediment of small grain size.This is consistent with ice samples retrieved from boreholes and trenches at this site.
Trench T7 at Site 2 (Figure 4b) shows black/brown 1961 tephra present as lapilli down to the ice table at 33 cm with some gray/brown ash from the same eruption.The 1875 tephra (Units C and D) are absent from the trenches at this site.The isopach maps from Carey et al. (2010) suggest that these layers may be present below the ice table we observed.The radargrams in Figure 3b show a reflector at the base of the ice in all cases.It is likely this reflector is generated by the massive ice in contact with these 1875 tephra units.We do not resolve multiple reflectors at depth, however.Unit C is likely thinner than Unit D at this site and is therefore not resolved by the radar antennas we used.

Site 3 Observations
At Site 3, both 1961 and 1875 tephra layers and buried ice are observed, including large hummocks of 1875 tephra and depressions infilled by 1961 tephra bounded by tension cracks (Figures 1e and 2c).GPR traverses were primarily taken along the banks of a large seasonal meltwater channel, which allowed easy access to vertical stratigraphy exposed through trenching.Trenches and boreholes from Site 3 are shown in Figures 4c-4f.
A trench (Figure 4c) dug at point T5 along the southwest traverse captures 1961 lapilli overlying 1875 tephra and massive ice.A zone of moisture extends from the top of the ice table to ∼20 cm within the 1875 ash layer.1875 tephra Units C and D are present in this layer.Figure 2c On the eastern bank of the channel, a trench revealed a transition from ∼25 cm coarse, mm-to-cm-size clasts of black/brown lapilli and gray/brown ash from 1961 to ∼15 cm of pore ice and 70 cm of massive ice.∼53 cm of an ice-cemented, pumice-supported coarse-to-fine ash matrix from 1875 extends to the floor of the channel (Units D and C).Another trench on the eastern bank ∼3 m of massive ice was buried beneath ∼40 cm of cm-sized 1875 tephra (Figure 4e).Pumice clasts in this image are from the 1875 Unit D (Carey et al., 2010).On the western channel bank within a large depression, ∼60 cm of ice buried by 36 cm of black/brown 1961 basaltic lapilli was uncovered on top of 1875 ash and pumice (Figure 4f).Both ice deposits had little entrained sediment and closed, mm-sized gas bubbles/pore space.
Figure 3c shows radargrams from D to D′ at Site 3, where up to 2.84 m of ice is buried exclusively by 62 cm of cm-to-dm-sized clasts of 1875 tephra and ash.The borehole retrieval shown in Figure 4e shows the transition from ice-cemented pumice clasts to massive ice.Moisture was minimal at this retrieval site.In contrast to the minimal scattering by the 1961 basaltic lapilli, the 1875 pumice layer shows strong scattering behavior in radargrams until it contacts the top of the ice table.A weaker basal interface is observed for the ice layer at both frequencies.The region between the ice table top and base is transparent in radargrams, suggesting fairly uniform ice with little entrained sediment or bedding at wavelength-scale, which was confirmed by both trenches exposing massive ice along the channel banks and a borehole taken close to D′ (see Figures 4d-4f).

Permittivity and Wave Velocity
Velocity measurements within the 1875 tephra were more challenging than those within the 1961 tephra.Scattering from 1875 pumice clasts approaching the wavelengths of the antennas generated many overlapping hyperbolic forms, especially at 400 and 900 MHz.Hyperbola fits were therefore scarce at Site 3, where the 1875 pumice is much more prevalent.We summarize all of the successful individual hyperbola fits in Table S1 and show examples of the fitting process in Figure S2 in Supporting Information S1.Average permittivity and velocities are primarily estimated from hyperbola fits for 1961 basaltic lapilli at each of the three sites and for each frequency.Hyperbola measurements from the 1875 tephra were much more uncertain.These permittivities are elevated compared to typical values of 2-4 for a dry, low-density volcanic tephra deposit (Campbell & Ulrichs, 1969).Averages are summarized for the hyperbola fitting method in Table 2.
Permittivity values estimated from trench and borehole measurements of tephra thickness are summarized in Table S2 in Supporting Information S1.We summarize averages of permittivity and corresponding velocity at all frequencies for all of the sites by capping the tephra deposit in Table 3.
Permittivity estimated from hyperbola fitting is elevated compared to the trench values.The hyperbola fitting method depends greatly on user identification of a true hyperbolic form and manual fitting of that shape (Giannakis et al., 2022).In particular, conventional hyperbola-fitting carries significant error when subsurface media are truly inhomogeneous (Giannakis et al., 2021).Fits to a false hyperbolic form yield unrealistic bulk permittivity values, which were likely retained in our sample leading to overestimates.This can be problematic even for automated hyperbola picking efforts where false and missed fits were shown to occur at 20% and 28%, respectively by Mertens et al. (2016).Furthermore, these hyperbola-fitting techniques can yield bulk velocities with variance and errors in the range of ±10% or more (Jol, 2008).We perform 10 repeat fits to 30% of all measured hyperbolas (N = 23) to estimate the variance and standard deviation of the velocity.We find our measurements to be self consistent based on our subset analysis with a low average variance and standard deviation of 2.03 × 10 4 and 0.01, respectively (see Table S4 for individual estimates of the subset).Fits to false hyperbolae are likely our greatest source of error.We report a relative error for permittivity estimated from Note.We report the standard deviation in the average permittivity calculated from the velocity.A relative error in the permittivity is reported compared to a laboratory measurement for a dry tephra.a See Campbell and Ulrichs (1969).
hyperbola fits in Table 2 and from trench measurements in Table 3 by comparing our estimates to a laboratory measurement where the maximum permittivity for a dry tephra was 4 (Campbell & Ulrichs, 1969).We primarily utilize permittivity estimates from trenches and boreholes in our analysis of losses and our ice thickness estimates since estimated relative error is consistently lower.

Calculated Losses to the Radar Signal
Based on our observations at each site, we assume a median permittivity of 9 for what are likely underlying basalt lava flows with some possible 1875 ash and/or tephra cover at Site 1 (Units C and D; Carey et al., 2010) and a value of 3 for a dry tephra to represent the 1875 ash and/or tephra underlying the ice layer along some traverses at Sites 2 and 3 (Campbell & Ulrichs, 1969) in order to perform a multilayer depth-correction prior to fitting.The permittivity assumption for Site 3 is likely an underestimate as we could not measure whether there was any water at the base of the ice layer or deeper below; substantial moisture would likely result in a higher bulk permittivity.Section 3.4 summarizes the permittivity assumptions for the ice and capping tephra layers.See Table S3 for a summary of all layer permittivities used for fitting.
Loss rate estimates result from portions of the corrected average trace (Figures 5a-5d, gray curves) that display a semilogarithmic decay with depth.In all cases, this portion of the curve was from the base of the ice table, or pumice layer in the case of the Vikrahraun lava flow, to the noise floor.The majority of the remaining sources of loss that are contributing to the estimated loss rates after the spreading correction are therefore below the ice table, or the less lossy pumice transitioning to the underlying 1961 Vikrahraun flow in the ice-free case.This is due to the fitting depth range of the corrected trace that spans from the sub-ice or sub-pumice layers to the noise floor.
The overlying tephra and ice contribute to the total loss, but likely to a lesser degree.We summarize the results of the attenuation for several representative ice-rich and ice-free radargrams in Figures 5a-5d at each of the antenna center frequencies.Average one-way losses are summarized for the various ice-rich and ice-free stratigraphy in Table 4. Figure S3 in Supporting Information S1 shows the average traces used for fitting in Figures 5a-5d.Table S3 summarizes individual loss rate calculations.
The ice-free 1961 Vikrahraun lava flow and ice-rich caldera sites share similar loss rates, particularly at 400 MHz.
There is some overlap between all calculated loss rates at 200 and 400 MHz in Figure 5e due to the spread in the values, likely the result of variability in scattering losses, soil moisture, or properties of the material beneath the ice between the different sites.Despite this overlap, there is an overall increase in one-way loss with increasing frequency, as observed by other investigations in permafrost environments (e.g., Boisson et al., 2011).Our loss rates are consistently lower than those measured by Boisson et al. (2011) at the same frequencies.Note.We report the standard deviation in the average permittivity calculated from the velocity.A relative error in the permittivity is reported compared to a laboratory measurement for a dry tephra.a See Campbell and Ulrichs (1969).b Only one traverse was taken at 400 MHz over purely 1875 tephra at Site 3 that had a coordinated borehole along that traverse (traverse 30, Figure 2c).Standard deviations cannot be estimated for this singular measurement.
We estimate fitting errors by examining the variance in the slope fit over the same depth range to subsets of 10 averaged traces at a time within the range of traces averaged for fitting.From this variance, we estimate the standard deviation for the one-way loss rates derived from the slopes of the fit, which are shown as error bars in Figure 5e.These are summarized in Table S3.Larger ranges in standard deviation result from poor fits and deviations from a semilogarithmic decay for each of the 10 traces sampled.There was no obvious correlation between the total number of traces and variance, indicating that higher variances result from variations in goodness-of-fit from each of the 10 trace samples.We found that the variance and standard deviation tend to increase with increasing frequency.The total number of loss rate estimates we were able to obtain was most often limited by a positive slope or low linear correlation coefficient, indicating a deviation from a semilogarithmic decay with depth.As a result, we obtained a limited number of loss estimates from Site 2 and radargrams with 1875 tephra burying ice (Site 3).

Bulk Radar Losses as an Indicator of Buried Ice
We tested whether total radar attenuation can be used as a diagnostic signature of buried ice to address the ambiguity that arises when interpreting radar data collected over stratigraphy sharing similar or overlapping dielectric properties.Askja was chosen for its substantial massive ice reservoirs buried by the dry, low-density, unconsolidated volcanic ash and tephra that directly overlap in permittivity with that of water ice (Campbell & Ulrichs, 1969;Johari, 1976;Matsuoka et al., 1997).Pure ice typically exhibits a lower loss than other geologic materials at microwave frequencies (Daniels, 2004;Pettinelli et al., 2015).The meters of ice in the caldera coupled with the tephra as a geophysical stand-in for regolith cover generated the best case scenario to test whether total attenuation yields a diagnostic signature of ice-rich versus ice-free stratigraphy.
Measured loss rates between ice-rich sites and the ice-free comparison site share much overlap in magnitude (Figure 5e), especially considering their associated formal errors (see Table S3).It is unclear whether the 1961 lava flow observation at 200 MHz in Figure 5e is anomalous.Perhaps the small grain sizes of the remobilized 1875 pumice and ash that blanket the lava flow relative to this wavelength combined with very dry soil create an extremely low loss environment in this case.In Figure S4 in Supporting Information S1, we compare our loss rates for ice-rich stratigraphy to literature values within the same range of frequencies.The tephra and ice at Askja may jointly result in a lower attenuation signature than the Fairbanks silt in Boisson et al. (2011).Boisson et al. (2011) note that the Fairbanks silt has a substantial clay content.Clays are a known radar attenuator due to adsorbed water and are problematic even in the most arid environments (e.g., Grimm et al., 2017).The lower total attenuation we estimate suggests that Askja offers a less attenuating subsurface than at Fairbanks, but it is challenging to isolate the signature of subsurface ice solely by evaluating trends in loss rates versus frequency.
Some of this ambiguity may be the result of the limited ice thicknesses encountered at our field locations and contributions of deeper volcanic materials to the measured loss rates.As discussed in Section 4.3, to achieve a reliable slope measurement, our one-way loss estimates require fits to the average trace over depth ranges that fall below the ice layer in all cases.Depending on the ice layer thickness, the combined loss of the tephra overburden and volcanic material that is below the ice (e.g., lava flows or older tephra layers) may contribute more greatly to the loss.Characterizing the scattering and absorption components of the total attenuation yield additional information about the subsurface and could aid in further isolating the signature of ice.Quantifying these components more robustly will be included in future analyses.If losses are instead plotted against average ice thicknesses by frequency, we find that the one-way losses decrease with increasing average ice thickness (Figure 6), suggesting that the signature of ice is likely to be more prominently expressed and identified as greater thicknesses are encountered.
Additional data points are needed to fully evaluate the relationship between one-way losses and average ice thicknesses.While all frequencies demonstrate this behavior, the 200 MHz data show a stronger correlation.The trends in Figure 6 are biased by the higher total number of loss rate fits that were able to be obtained at 200 MHz.

Other Sources of Loss
The total losses we estimated are comprised of contributions from both absorption and scattering.Without a full characterization of these effects, it is difficult to determine if losses are simply driven by the different scattering regimes created by the various tephra clast sizes (mm-scale 1961 lapilli vs. cm-scale 1875 pumice) and material below the ice layer (such as lava flows or older tephra), or absorption, or a combination of both.In Figure 5, for example, the larger cm-to-dm-scale 1875 pumice clasts overlying massive ice exhibit a lower loss at 200 MHz than the observations over the 1961 lapilli at the same frequency.We observed that scattering became more visually prevalent in the 400 and 900 MHz radargrams in the 1875 pumice deposit, suggesting that scattering is dominant at these frequencies and to a lesser extent at 200 MHz.In the 1961 basaltic tephra deposit, absorption may be more dominant at 200 MHz given the sub-wavelength clast sizes (lapilli) associated with that deposit.Additionally, traverses with 1875 pumice coverage have thicker ice layers on average.Given the correlation between lower loss rates and average ice thickness (Figure 6), the lower loss rates observed for the 1875 pumice traverses could be the result of a thicker ice layer at depth for those traces.This highlights the need for further independent measurements and modeling of contributions from absorption and scattering, respectively, in order to distinguish an ice-rich from an ice-free subsurface.
Figure 6.One-way losses decrease with increasing average ice thickness along sections of ground-penetrating radar traverses comprising the average trace at each site.This trend is likely biased by the higher number of loss rate fits obtained at

Insights on Planetary GPR Deployment
Askja presents a testbed to conduct controlled geophysical experiments for future operations, instrument deployment, and post-collection analyses that can be used to identify resources of interest for real-time decisionmaking during surface missions.The unvegetated terrain within the caldera and surrounding regions share a striking similarity to the surfaces of the Moon and Mars.Unsurprisingly, the ambiguity associated with GPR returns from the subsurface there necessitated additional context such as boreholes or trenching in order to confirm our interpretation of stratigraphy.This will no doubt be a challenge on the surfaces of other terrestrial bodies for both robotic and human assets.Here, we discuss lessons learned from Askja and implications for GPR field deployment and analysis of radar attenuation for future ISRU analog investigations and more broadly cryogeophysics in volcanic terrains.
The frequency range used to map the tephra and ice in the caldera was successful at resolving the entire vertical extent of both the tephra and massive ice, where present, at all sites.All three center frequencies performed exceptionally well when taking data over the mm-sized 1961 basaltic lapilli.200 MHz was particularly advantageous for overcoming the effects of scattering from the cm-to-dm-sized clasts of 1875 pumice, especially at Site 3. 400 MHz provided a better balance of vertical resolution and penetration depth, particularly at Sites 1 and 2. Surveying at the 900 MHz center frequency was challenging at all sites except for Site 1, where a combination of conditions including smoother, more compact tephra (enabling constant ground coupling) and small clast sizes allowed for the collection of high-quality radargrams.The finer vertical resolution at this frequency also enabled greatly improved identification of bedding and internal layering in the tephra.For any resource campaigns conducted on the surface of the Moon or Mars where higher amounts of scattering in the regolith may be a concern, lower frequencies will likely be more favorable, unless the characterization of very near-surface ice is desired, such as what was encountered at the high latitudes of Mars by the Phoenix Lander (Mellon et al., 2009).
Our GPR surveys across the three caldera floor sites did not resolve vertical transitions between pore ice and massive ice.Our trench and borehole observations showed thin (generally <10 cm thick) zones of pore-ice within the tephra before transitioning to massive ice at depth.These zones in pore ice and transitions in ice concentration are at or below the limit of vertical resolution for the GPR frequencies used here, making them difficult to identify in radargrams.To test GPR's ability to characterize changes in ice concentration for future planetary science and ISRU applications, more suitable analog sites should be sought.This will be important for future interpretation of Mars ground data and especially for the Moon, since ice in the form of pore-filling zones or small grain sizes may be most prevalent, at least in the near-surface (Siegler et al., 2015).
The Askja caldera tephra deposits are not precisely compositionally analogous to lunar or martian regolith but are of a similar density and permittivity when dry (see Carrier et al., 1991;Lai et al., 2019;Olhoeft & Strangway, 1975).Moisture, although relatively minimal at Askja at the time of surveying, will be an issue for GPR investigations at any Earth analog site.This is the key difference between arid, present-day Mars, the Moon, which is an airless body, and any Earth analog site.Although the effects of this moisture on the attenuation are likely minimal at Askja, this thin zone of elevated moisture at the top of the ice table likely enhanced the dielectric contrast between the ice and overlying tephra, which are otherwise close or overlap in permittivity (e.g., Boisson et al., 2011;Campbell & Ulrichs, 1969).On Mars, the contact of layers of similar permittivity may not produce a reflector at all in such a dry environment.In particular, permittivity cannot be used alone to identify buried water ice on Mars, where it can easily be mistaken for a low-density volcanic tephra or ash.An excellent example of this ambiguity is the Medusae Fossae Formation, where an ashfall deposit or water ice were considered equally likely sources of observed permittivity and layers from orbital radar data.A recent investigation utilizing a combinatination of permittivity, radar loss estimates, and modeling has revealed this deposit is likely to be layers of dust and water ice (Watters et al., 2024).This highlights the challenge of uniquely identifying water ice using permittivity alone.Advanced GPR measurement and analysis techniques using the phase characteristics and polarimetric response of returned signals not explored by our investigation could also be assessed for their utility in identifying and characterizing subsurface ice.
Moisture at our field sites likely contributed to the bright reflectors observed in Figure 3 at the tephra-ice table contact.Characterization of moisture at any analog site is necessary in order to isolate its contribution to radar attenuation.In order to directly compare results from future terrestrial analog field GPR investigations with the arid martian or lunar environments, soil moisture probe measurements with depth or coordinated electrical resistivity measurements to quantify contributions to the total attenuation from meltwater or moisture should be included in the field plan.We attempted to take soil moisture measurements with a probe in 2021 and observed general trends of increasing moisture with depth toward the ice table, but they were challenging to make accurately due to the coarse-grained nature of the tephra.Due to the potential inaccuracies this introduces, we do not quantify how it contributes to attenuation.A complementary method to overcome this challenge could be the measurement of gravimetric water content.Despite this, our loss rates are lower than those observed in the Fairbanks silt, which is known to have a significant clay content that is attenuating to radar (Boisson et al., 2011;Grimm et al., 2017).The effects of moisture may be further minimized by conducting field investigations at times of year when the ground is still frozen though safe access may be more challenging.Other options include investigations in extreme environments such as the high Andes or the Antarctic Dry Valleys.Despite these limitations in environmental conditions, the Askja caldera deposits offer variation in clast sizes, intra-grain and intergrain porosities, varying ice thicknesses, and ice burial by a low-density overburden.This diversity makes Askja not only an ideal testbed for ISRU applications but also an important site for studying volcano-ice interactions and Icelandic permafrost reservoirs.

Conclusions
We conducted the first GPR survey of tephra sourced from the eruptions of Askja in 1875 and 1961 that buried and preserved layers of massive ice.The overlapping dielectric properties of the Askja tephra deposits and metersthick massive ice layers make this site an ideal location to test GPR data analysis techniques for future ISRU objectives for the Moon and Mars.Determining detection thresholds for water ice deposit thicknesses and concentrations is important for future interpretation of GPR observations at the Moon and Mars where concentrations are expected to be lower and ice-regolith mixtures could produce a nonunique signature.We estimated and compared total attenuation rates between ice-rich and ice-free stratigraphy to determine if massive ice at depth produces a diagnostic loss signature at multiple GPR center frequencies.
GPR successfully delineated the vertical and horizontal extent of shallow ice deposits at Askja.Loss rates across all frequencies and at all sites with different tephra cover and ice thicknesses were consistently low.Loss rates increase with increasing frequency and decrease with increasing ice thickness.Thicker layers of ice that comprise a larger fraction of the overall attenuation path length of the radar signal may yield improved detectability using this method.
While attenuation may not be solely diagnostic of the signature of water ice, when combined with analyses of reflectors and hyperbola fitting, a unique result can be obtained.Pore ice, ice and regolith mixtures, and generally low concentrations of ice will be challenging to uniquely identify using traditional approaches to GPR data analysis.More data points at additional frequencies would provide more confidence in these trends.Additional modeling and measurements of contributions to the total attenuation from scattering and absorption will be necessary to further isolate the signature of massive ice.

Figure 1 .
Figure 1.Askja is a central volcano located in Northern Iceland (panel a inset).(a) Deposits of buff-colored rhyolitic pumice and black/brown basaltic tephra from eruptions in 1875 and 1961, respectively, cover the northeastern part of the caldera.These deposits preserve ice at depth.Ground-penetrating radar (GPR) surveys were conducted within the caldera and at the toe of the 1961 Vikrahraun lava flow (red boxes).(b) Within the caldera, GPR surveys were concentrated in two main regions: within the 1961 tephra deposit to the north (panel c) and a zone where 1961 tephra transitions to 1875 pumice (panels d, e).(c, d, e) Sites 1, 2, and 3, respectively, where GPR data were collected at 200, 400, and 900 MHz in summer 2019 and 2021.GPR traverses at all center frequencies taken in 2019 and 2021 are collectively shown in yellow.Panel (a) basemap is Landsat-8 pansharpened image LC8_217015_20140906 (Vermote et al., 2016).Panel (b) basemap obtained from the Esri ArcGIS™ software.

Figure 2 .
Figure 2. We mapped massive ice buried by tephra using ground-penetrating radar (GPR) at center frequencies of 200, 400, and 900 MHz.Panels a, b, and c are ice thickness maps summarized for the 400 MHz observations.Massive ice was observed to be thickest at Site 2, panel (b), reaching depths of 4.4 m.Ice thicknesses <10 cm are reported as zero on these maps.Radargrams collected at each center frequency at each of the three sites are summarized in Figure 3. Trench and borehole locations are marked in red.(a) Site 1 massive ice deposits are thinner in the west and thicker toward the hiking trail in the east, reaching a maximum depth of 2.8 m.The tephra transitions from more loosely packed ripples or aeolian bedforms to a smoother, compact section after point A. Internal layering was observed within the less compact tephra.(b) Site 2 has more ablation features than Site 1 (e.g., tension cracks and depressions).Massive ice is likely resting on an older ash or tephra deposited prior to 1961.Massive ice deposits are also thickest at this site, reaching a maximum of 4.4 m moving south along the hiking trail.(c) Massive ice at Site 3 tends to be thicker, on average, beneath the 1875 pumice.Obvious thermokarsts are scattered throughout this region.A meltwater channel generated by seasonal melting of snow cuts through GPR traverses taken on either side.Trenches taken on the east and west banks of this channel revealed ice buried by tephra from the 1875 and 1961 eruptions.Images of the interplay of these deposits and massive ice along with labeled trenches and boreholes from this figure are shown in Figure 4.

Figure 3 .
Figure 3. (a) Site 1 radargrams correspond to panel (a) of Figure2.A section of a 206 m traverse is shown from A to A′ where up to ∼2 m of massive ice occupies the troughs of buried lava flow festoons.The top of the ice table is fairly constant, measured from trenching to be at 52 cm depth beneath the 1961 tephra.From the ice table position, an average permittivity of 2.91 is applied to place the ice table at its observed depth in the images.(b) Site 2 radargrams correspond to panel (b) of Figure2.Section B to B' is a 900 MHz observation with relatively thin ice (∼50 cm) preserved by ∼33 cm of 1961 tephra taken in 2021.This radargram was depth-corrected using a derived tephra permittivity of 7. Section C to C′ is repeat 400 and 200 MHz observations taken in 2019 where up to 65 cm of massive ice is preserved by ∼50 cm of 1961 tephra.These radargrams were depth-corrected using a permittivity of 2.59.(c) Site 3 radargrams correspond to panel (c) of Figure2.Radargrams from D to D′ were taken upslope.The cm-to-dm-sized clasts of pumice from 1875 are visibly more scattering than the sub-wavelength 1961 tephra.Massive ice reaches a maximum thickness for Site 3 along this section at 2.8 m.The ice table is estimated to be at ∼62 cm from a borehole 29 m along-traverse.These radargrams were depth-corrected using an average tephra permittivity of 4.6.All radargrams were processed in the GSSI RADAN 7 software using a vertical infinite impulse response filter, a five-point exponential gain, and corrected for variations in surface topography.

Figure 4 .
Figure 4. Trenches and boreholes were taken along various ground-penetrating radar (GPR) traverses at each of the three sites in the caldera and on the 1961 Vikrahraun lava flow.Trench and borehole numbers correspond to labeled red points in Figure2.Sites 1 (a) and 2 (b) are blanketed by ∼30-50 cm of black/brown basaltic lapilli mixed with some gray/brown ash from the 1961 eruption overlying massive ice.These trenches were dug to the top of the ice table in each case.1875 tephra (Units C and D) are likely present below the ice table overlying other older basalt lava flows (e.g.,Carey et al., 2010) but were not observed in the trenches.Site 3 trenches in panels (c) and (d) show the interplay between black/brown 1961 tephra, massive ice, and the 1875 tephra deposits.In panel (c), ∼15 cm black/brown lapilli and coarse ash from 1961 give way to 15 cm of moist brown ash with mm-to-cm-sized blocks of buff colored 1875 pumice.The ice table rests at 30 cm depth.T2 in panel (d) was dug on a slope of the eastern meltwater channel bank in Figure2c, making observed depths larger than vertical.∼25 cm coarse, mm-to-cm-size clasts of black/brown 1961 tephra and gray/brown ash transitions to ∼15 cm of pore ice and 70 cm of massive ice.∼53 cm of an ice-cemented, pumice-supported coarse-to-fine ash matrix from 1875 extends to the floor of the channel, which is comprised of a 50/50 mixture of pumice clasts and ice from 1875.(e) Borehole B7 (Figure2c) was taken at 29 m along radargrams in Figure3c.This retrieval (arranged in discontinuous but depth-ordered sections) shows the transition from ice-cemented 1875 pumice and ash to massive ice from bottom to top in the image.The ice table was measured to be at 62 cm depth.The trench was taken along the eastern meltwater channel bank near these GPR traverses.∼40 cm of cm-sized 1875 pumice with ash is burying ∼3 m of massive ice at Site 3. (f) This trench was dug along the western bank of the meltwater channel at Site 3 where 36 cm of 1961 basaltic lapilli was covering ∼60 cm of massive ice that was resting on moist, brown 1875 ash deposits with some mm-to-cm-sized clasts of buff colored 1875 pumice.(g) GPR traverses across the toe of the 1875 pumice-mantled 1961 Vikrahraun lava flow found ∼1 m of pumice rests atop the basalt flow.1875 tephra was finer-grained at this site (mm-sized grains and fine-grained ash).Stratigraphy is ice-free.

Figure 5 .
Figure 5.We fit sections of semilogarithmic decay of trace averages at various sites and frequencies across the Askja caldera (a-c) and from the 1961 Vikrahraun basalt lava flow mantled with 1875 pumice (d).Figure S3 in Supporting Information S1 shows the traces averaged in this figure on their respective radargrams.The average trace (solid black line) is corrected for geometric spreading effects (gray solid line).Sections of semilogarithmic decay with depth are fitted with a linear regression (blue dashed line).Losses are comparable between the Site 3 and 1961 Lava Flow observations at the same frequency (Panels c and d).Thicknesses of the shaded regions correspond to the estimated thickness of the deposit for those averaged traces.The thickness of the deposit(s) underlying the ice is unknown in all cases.One-way loss in dB/m calculated from the slope of the fit (blue-dashed line), Q, and the loss tangent (tan δ) are shown.(e) One-way attenuation at 200, 400, and 900 MHz increases with increasing frequency.Points at each center frequency are spread out for clarity.There is a significant spread in the losses associated with each deposit at each frequency.We report the slopes of successful fits along with other parameters in TableS3.
Figure 5.We fit sections of semilogarithmic decay of trace averages at various sites and frequencies across the Askja caldera (a-c) and from the 1961 Vikrahraun basalt lava flow mantled with 1875 pumice (d).Figure S3 in Supporting Information S1 shows the traces averaged in this figure on their respective radargrams.The average trace (solid black line) is corrected for geometric spreading effects (gray solid line).Sections of semilogarithmic decay with depth are fitted with a linear regression (blue dashed line).Losses are comparable between the Site 3 and 1961 Lava Flow observations at the same frequency (Panels c and d).Thicknesses of the shaded regions correspond to the estimated thickness of the deposit for those averaged traces.The thickness of the deposit(s) underlying the ice is unknown in all cases.One-way loss in dB/m calculated from the slope of the fit (blue-dashed line), Q, and the loss tangent (tan δ) are shown.(e) One-way attenuation at 200, 400, and 900 MHz increases with increasing frequency.Points at each center frequency are spread out for clarity.There is a significant spread in the losses associated with each deposit at each frequency.We report the slopes of successful fits along with other parameters in TableS3.

Table 1
Estimates of Average and Maximum Massive Ice Thickness, Average and Maximum Tephra Overburden Thickness, and Their Standard Deviation From 400 MHz Ground-Penetrating Radar Observations for the Three Caldera Sites SHOEMAKER ET AL.

Table 2
Average Permittivity and Velocity From Hypberbola Fits for the Three Caldera Sites

Table 3
Summary of Average Permittivity and Velocity for Each Tephra Deposit Derived From Trenches and Boreholes

Table 4
Summary of Average One-Way Attenuation for Ice-Rich Caldera Stratigraphy and the Ice-Free, 1875 Pumice-Mantled 1961 Vikrahraun Lava FlowNote.Standard deviation is reported for attenuation estimated for each frequency and inferred stratigraphy.
SHOEMAKER ET AL.