Age-Depth Stratigraphy of Pine Island Glacier Inferred From Airborne Radar and Ice-Core Chronology

However, much less is known about the interior ice-sheet history of this region during the Holocene. Cosmogenic Abstract Understanding the contribution of the West Antarctic Ice Sheet (WAIS) to past and future sea level has been a major scientific priority over the last three decades. In recent years, observed thinning and ice-flow acceleration of the marine-based Pine Island Glacier has highlighted that understanding dynamic changes is critical to predicting the long-term stability of the WAIS. However, relatively little is known about the evolution of the catchment during the Holocene. Internal reflecting horizons (IRHs) provide a cumulative record of accumulation, basal melt, and ice dynamics that, if dated, can be used to constrain ice-flow models. Here, we use airborne radars to trace four spatially extensive IRHs deposited in the late Quaternary across the Pine Island Glacier catchment. We use the WAIS Divide ice-core chronology to assign ages to three IRHs: 4.72 ± 0.28, 6.94 ± 0.31, and 16.50 ± 0.79 ka. We use a 1-D model, constrained by observational and modeled accumulation rates, to produce an independent validation of our ice-core-derived ages and provide an age estimate for our shallowest IRH (2.31–2.92 ka). We find that our upper three IRHs correspond to three large peaks in sulfate concentrations in the WAIS Divide ice-core record and hypothesize that the origin of these spatially extensive IRHs is from past volcanic activity. The clear correspondence between our IRHs and the ones previously identified over the Weddell Sea Sector, altogether representing ∼ 20% of the WAIS, indicates that a unique set of stratigraphic markers spanning the Holocene exists over a large part of West Antarctica.

nuclide studies on isolated nunataks across the ASE suggest significant ice thinning occurred during the early-to mid-Holocene in the central ASE (Johnson et al., 2017(Johnson et al., , 2020Lindow et al., 2014), with thinning complete by the mid-Holocene in the eastern ASE near PIG (Johnson et al., 2008. More recent evidence, based on sediment cores, ice-penetrating radar, and ice-sheet modeling, showed possible retreat and re-advance of the WAIS grounding line over millennial timescales during the Holocene (Kingslake et al., 2018), although evidence of such behavior is not available in the ASE region.
Internal reflecting horizons (IRHs), as observed by ice-penetrating radars, provide a powerful and complementary resource to point-based geochronological measurements. Excluding basal ice and erosional surfaces, the majority of specular, continuous IRHs are isochronous (Whillans, 1976); many can be traced for several hundreds of kilometers and provide a record of accumulation rates and patterns, convolved with key information on past ice-dynamical processes (Bingham & Siegert, 2007;Eisen et al., 2005Eisen et al., , 2008Siegert et al., 1998). IRHs can thus serve as a valuable resource for constraining past changes in surface mass balance (SMB) and ice-flow velocities (e.g., Rotschky et al., 2004), and, where they can be dated, can be incorporated into ice-flow models, as previously shown for Greenland MacGregor et al., 2016) and Antarctica (Cavitte et al., 2018;Koutnik et al., 2016;Leysinger Vieli et al., 2011;Waddington et al., 2007).
Despite the large spatial coverage of radar data across Antarctica, information on dated IRHs is limited over much of the WAIS. This is partly due to the restricted availability of deep ice cores, the multitude of radar-system families operating at varying frequencies and using different post-processing methods to generate the radar data, and the challenge in tracing deep continuous IRHs, particularly through areas of high strain rate (i.e., at the onset of fast-flowing tributaries). Nonetheless, previous studies over the WAIS have used IRHs for the direct purpose of linking major deep ice cores Neumann et al., 2008), while others have used a wider, catchment-scale approach to constrain information on past accumulation rates and ice-flow reconfiguration. Such studies have ranged across the central WAIS (Jacobel & Welch, 2005;Muldoon, 2018;Siegert & Payne, 2004) or focused on specific sub-regions, for example, Siple Dome (Jacobel et al., 1996), Kamb Ice Stream (Catania et al., 2006;Holschuh et al., 2018) and Thwaites Glacier .
Over PIG, Karlsson et al. (2014) identified two IRHs spanning much of the slow-flowing parts of the catchment, which they roughly dated to 5.3-6.2 and 8.6-13.4 ka. More recently, Ashmore et al. (2020) recovered three IRHs ranging across Institute and Möller Ice Streams and crossing the Institute/PIG divide which they broadly dated to 1.9-3.2, 3.5-6.0, and 4.6-8.1 ka. They demonstrated a correspondence between their IRH package and the IRHs previously identified by Karlsson et al. (2014) and Siegert et al. (2005), suggesting that a spatially extensive network of IRHs may span much of the WAIS.
Here, we build on previous studies to present a spatially extensive, dated-radiostratigraphy of PIG. We use ice-penetrating radar data collected from two airborne platforms to trace four IRHs throughout PIG. We use a published ice-core chronology as well as a steady-state vertical-strain model to date these IRHs, and show that they span much of the Late Pleistocene and Holocene. We first discuss the specifications of the radar systems and their respective uncertainties, and then describe the methods used to assign ages to each of our four IRHs. We present the dated age-depth stratigraphy of the catchment and make inferences for the rest of WAIS by comparing our recent findings to other age-depth studies. Finally, we investigate the link between sulfate activity in the WAIS Divide ice-core record and the depth of our upper three IRHs, and discuss the potential to recover records of older (i.e., pre-LGM) ice in the region using currently available radar data sets.

Data
The principal data used in this study were acquired during two large-scale airborne radar surveys of West Antarctica.
The first of these was acquired over the 2004-05 austral season, when PIG's 175,000 km 2 catchment was surveyed extensively using the British Antarctic Survey's Polarimetric Airborne Survey INstrument (PASIN) BODART ET AL. 10.1029/2020JF005927 2 of 20 system . This survey, hereafter termed "PIG-PASIN", acquired ∼35,000 line-km of airborne radar data across the region (Figure 1). Data were collected with two interleaved radar modes. The first was a deep-sounding, 150 MHz center-frequency, 4-µs, 10-MHz chirp mode, which has been used previously to identify and trace the bed  and some IRHs (Karlsson et al., 2009(Karlsson et al., , 2014. The second was a 150 MHz, 0.1-µs pulse mode designed to image shallow IRHs but from which we are also able to recover IRHs deeper (∼2 km, see Figure 2a) in the ice column. Over much of the surveyed region, flight lines form 30 km spaced grids that contain multiple crossovers, ensuring consistency when tracing IRHs across neighboring lines (Figure 1). Following techniques outlined in Ashmore et al. (2020), here, we used both modes of PASIN interchangeably during our IRH-tracing procedures (see Section 2.2). For the purposes of linking our stratigraphy further across the WAIS, we also refer to further PASIN-acquired data from a survey of Institute and Möller Ice Streams undertaken in 2010-11 (hereafter "IMAFI-PASIN"), BODART ET AL.

10.1029/2020JF005927
3 of 20 Figure 1. Map of study area with the data sets and key locations mentioned in this study. The inset in top left corner shows the region of interest (red box). Airborne survey lines included in this study: PIG-PASIN (gray), OIB-MCoRDS2 (yellow), IMAFI-PASIN transects flown over Institute Ice Stream (IIS) and intersecting the PIG catchment (orange), SPRI/NSF/TUD line (brown), overlaid on top of ice flow velocities from Rignot et al. (2017) and MODIS Mosaic of Antarctica (Scambos et al., 2007). Also included is the long, ITASE GPR-transect (dashed red) through which the 17.5 ± 0.5 ka layer from Jacobel and Welch (2005) was traced. The numbers shown over PIG's trunk represent the eight fast-flowing tributaries (1)(2)(3)(4)(5)(6)(7)9) mentioned in this study. The WAIS Divide (WD2014) and Byrd ice cores are represented by the two black triangles, and the black arrows represent the three intersections between the SPRI/NSF/TUDtraced IRHs and this study. The two red circles show the two sites (Sites A and B) where the 1-D age-depth model was used. The AA-AB segment (magenta) shows a subset of the control line where IRHs were first identified over PIG-PASIN (see Figure 2). The Western Divide is shown as WD on the map. The ICESat IMBIE basins of Pine Island Glacier (PIG) and Thwaites Glacier (TG) (Zwally et al., 2012) are annotated on the map and delimited by the blue outline lines. which provided tie-lines connecting PIG with its neighboring basins (Figure 1; see Ashmore et al., 2020, and references therein, for further details).
The second survey was conducted in 2016 and 2018 by NASA's Operation IceBridge (OIB) mission, and yielded ∼3,000 line-km of airborne radar data over PIG, Institute and Möller Ice Streams, and Thwaites Glacier ( Figure 1). The system deployed by the Center for Remote Sensing of Ice Sheets (CReSIS) was the Multichannel Coherent Radar Depth Sounder 2 (MCoRDS2) with a 190 MHz center frequency and 50 MHz bandwidth. We used the CReSIS L1B standard products, produced with pulse compression, focused-SAR processing, and along-track motion compensation. More information on the radar system and processing is given by CReSIS (2016). Critically for this study, one of the OIB flight tracks over PIG also flew over the WAIS Divide Ice Core (79.48°S, 112.11°W; hereafter referred to as WD2014) (Figure 1), making it possible to assign relatively unambiguous dates to the traced IRHs.
More details on each of the radar systems are provided in Table 1. For the purposes of increasing IRH traceability on the PIG-PASIN data, we quadratically detrended each radar trace, normalized each pixel in a moving vertical window, and then applied a 10-trace horizontal average to reduce incoherent noise (after Ashmore et al., 2020). For both the PIG-PASIN and the OIB-MCoRDS2 data, we removed the air-to-ice two-way travel time and shifted the surface elevation to time zero, before exporting the data to standard 2-D SEG-Y format for data interpretation.

IRH-Tracing Workflow
We conducted all IRH-tracing in the Schlumberger Petrel ® 3-D seismic software using a semi-automated tracing algorithm that uses an adjustable window to track the local maxima of received reflected power between traces.
We initiated our workflow on the PIG-PASIN data set as it is the most spatially extensive survey of the PIG catchment. From a "control line" crossing the ice divides between PIG, Thwaites Glacier, and Institute Ice Stream (Figure 1), in which clearly visible englacial stratigraphy is ubiquitous in both chirp-and pulsemode data, we identified four prominent IRHs that we term R1-4 ( Figure 2). The upper three IRHs (R1-3) were chosen on the basis of high spatial continuity, high signal-to-noise ratio (SNR), and as being analogous to "IRH packages" traced over part of PIG by Karlsson et al. (2014) and through IMAFI-PASIN radar profiles by Ashmore et al. (2020). All four IRHs occur in the middle part of the ice column where IRHs likely result from contrasts in acidity from past volcanic eruptions (Gow & Williamson, 1971;Millar, 1981Millar, , 1982, rather than the result of density variations occurring primarily at the near-surface (Clough, 1977;Gow, 1970;Moore, 1988) or orientation of anisotropic material due to ice foliation in the basal zone (Fujita et al., 1999;Harrison, 1973); and thus can be assumed to be isochronous (Siegert et al., 1998;Whillans, 1976).
Expanding out from the control line, we progressively traced and mapped IRHs across the catchment using IRH intersections at each crossover as calibration points. This ensured reliability in our reflection tracing as the software is capable of detecting intersecting IRHs at the crossover with orthogonal radar lines. Since our tracing strategy was based on reflector echo strength and continuity, the reflection tracing was terminated when it was no longer possible to distinguish visually between adjacent reflections, either as a result of similar brightness levels or a loss in continuity. This was particularly common in areas of steep bed topography BODART ET AL. Note: for the PASIN system, we provide values for both the chirp-and pulse-acquisition mode in the bandwidth/ pulse width column, as well as in the vertical resolution column. The vertical resolution of the chirped systems was calculated as per CReSIS (2016) using a scaling factor "k" which accounts for resolution degradation due to receiver characteristics and processing (see Equation S1 ).

Table 1 Characteristics and Resolution of the Two Airborne Radar Systems Used in This Study
causing IRHs to dip significantly, or where enhanced ice-flow speeds disrupted IRH continuity, notably into the main flow features of PIG's northern catchment. In some places, IRHs faded without such clear topography/flow-induced reasons, likely due to the attenuation of the radar signal with depth or the type of processing used (Holschuh et al., 2014). In some locations more distant from the upper PIG catchment (i.e., westward of tributary 6; Figure 1), extensive englacial layering was visible in radar profiles but, due to a dearth of connecting lines and crossovers, we could not, with confidence, identify R1-4.
When tracing between crossovers, we relied upon the distinctiveness of our IRHs. At the vertical resolution of PASIN, R1 and R2 manifest as single-amplitude peaks, with R2 representing a particularly bright reflector widely visible across our radar data ( Figure 2, Figure S1). R3 consists of the shallowest of a series of closely spaced bright horizons, often manifested as a couplet (zoomed inset in Figure 2, Figure S1), and previously identified by Karlsson et al. (2014; their "Layer 2") and Ashmore et al. (2020;their "H3"). The lowermost IRH, R4, forms the upper part of a band of bright reflectors visible at the intersection with the 17.5 ± 0.5 ka layer widely imaged on radar data from the International Trans-Antarctic Scientific Expedition (ITASE) connecting the PIG catchment with the Byrd Ice Core chronology (Hammer et al., 1997;Jacobel & Welch, 2005) (Figures 1, 2, and S1).
Once R1-4 were traced through the PIG-PASIN survey, we looked for the same IRHs on the OIB-MCoRDS2 data using available crossovers between each survey (Figures 1 and 3). We found R2-3 to be equally distinguishable in OIB-MCoRDS2 profiles, with R2 representing a particularly bright reflector similar to that on PIG-PASIN, whilst R3 also formed the shallower part of an easily distinguishable couplet. We did not recover R1 independently on the OIB-MCoRDS2 profile crossing the WAIS Divide Ice Core and used intersections with PIG-PASIN to trace it across to the Institute Ice Stream catchment. Similarly, we used several intersections with the 17.5 ± 0.5 ka layer from Jacobel and Welch (2005) in and around the WD2014 site to recover R4 in the OIB-MCoRDS2 data (Figures 1 and 3).
It is worth noting that the OIB-MCoRDS2 data were acquired 12-14 years later than the PIG-PASIN survey, and so the same IRHs will, in principle, lie slightly lower in the ice column. However, considering a present-day mean accumulation rate of ∼0.30-0.35 m a −1 (meters of ice equivalent per year) at the intersection between the two surveys, the maximum change in IRH depth is < 5 m. This is well within the bounds of the total depth uncertainty calculated for each radar system (see Section 2.3) and does not affect the pattern BODART ET AL.
10.1029/2020JF005927 5 of 20 of englacial layering or the identification of our IRHs across the different surveys. Crossover analysis at key intersections on the airborne radar data showed that the mean depth difference for R1-4 falls within the uncertainty range of all surveys ( Figure S2, Tables S1 and S2) (see Section 2.3). At 10 intersections on PIG-PASIN, the mean depth difference for R1-4 is < 6 m. Similarly, the mean depth difference for R2-3 at 11 intersections between PIG-PASIN and OIB-MCoRDS2 is 14 and 29 m, respectively, and < 18 m at 5 intersections between R4 on OIB-MCoRDS2 and the 17.5 ± 0.5 ka from Jacobel and Welch (2005) ( Figure S2, Table S2).
With our objective being to produce an age-depth radiostratigraphy across PIG, we converted all IRHs traced above in the time domain (t IRH ) to depth (d IRH ) using is the speed of electromagnetic waves through ice (cf. Fujita et al., 2000) and Z f = 10 m is a spatially invariant firn correction, appropriate for West Antarctica (Ashmore et al., 2020). All our depth measurements are given in meters below the surface. We then calculated IRH depth as a function of ice thickness using the ice-thickness measurement from each respective radar mission, and complemented these with ice-thickness measurements from BedMachine (Morlighem et al., 2020) in places where the radar did not sound the bed.

Catchment-Wide Depth Uncertainties
To assess the accuracy of our IRH depths at the catchment scale, we consider the uncertainties associated with the imaging of IRHs with ice-penetrating-radar. These uncertainties primarily depend on three factors: variations in the speed of electromagnetic (EM) waves through the ice, the firn-density correction, and the radar system's range precision (Cavitte et al., 2016) (Text S1).
The maximum uncertainty arising from selecting an EM value ranging between 168 and 169.5 m μs −1 is 16 m on the maximum depth of the deepest reflection on PIG-PASIN and 14 m on OIB-MCoRDS2. The uncertainty associated with the firn correction is ±3 m, owing to minor variations in firn densification across the catchment (Ashmore et al., 2020) (Text S1). The precision of IRH depth estimates also depends on the range accuracy, σ(r * ), of the radar system, which refers to how accurately changes can be located in 3-D space (Cavitte et al., 2016;King, 2020). This is a combination of the SNR of each IRH and the range resolution, Δr, of the radar system, which is mainly a function of sampling frequency, bandwidth, source wavelets, and the type of post-processing applied. The range resolution for each system, from coarser to finer is: PASIN chirp (12.89 m), PASIN pulse (8.42 m), and MCoRDS2 (2.58 m) ( Table 1, Text S1).
We undertook an empirical error analysis to calculate the maximum uncertainty associated with the deepest IRH by calculating the root-mean-square error of the depth uncertainties from EM wave through the ice, the firn correction, and the radar range accuracy. We obtained a combined maximum uncertainty of ±17 m and attached this uncertainty to all IRHs traced on the PIG-PASIN data (Text S1). Similarly, we estimated a combined maximum uncertainty of ±14 m on the OIB-MCoRDS2 data (Text S1). Given that this uncertainty represents the maximum uncertainty on the deepest IRH over our entire data set, we also calculate IRH-specific uncertainties at the ice-core site (see Section 2.4.1).

Age-Depth Attribution
To estimate the absolute age of our IRHs, we employ two primary dating methods: we use (a) the WAIS Divide ice-core chronology to provide a direct age to our three deepest IRHs, namely R2-4; and (b) the Dansgaard-Johnsen 1-D model to independently compare the ages calculated at the ice core and to provide an approximate age range to our shallowest IRH, R1. Once dated, we also compared the ages and depths of R1-3 with dated IRHs traced across PIG (Karlsson et al., 2014;Siegert & Payne, 2004) and Institute and Möller Ice Streams (Ashmore et al., 2020); as well as the age and depth of R4 with the 17.5 ± 0.5 ka layer dated using the Byrd ice-core chronology (Hammer et al., 1997)

Connection to the WAIS Divide Ice-Core Chronology
We used the 2016 OIB-MCoRDS2 data linking central PIG to the WD2014 site to date IRHs across PIG relative to the ice-core chronology, where annual-layer counting goes back to the last ∼31 ka BP (Buizert et al., 2015;Sigl et al., 2016). We take the recorded depth at the ice core which most-closely matches our IRH depth at WD2014, and calculate the upper and lower age bounds using the radar depth and ice-core BODART ET AL.
10.1029/2020JF005927 7 of 20 uncertainties. Following MacGregor et al. (2015), the age uncertainty (Δa comb ) associated with each IRH is the root-mean-square combination of the age uncertainty associated with the unweighted mean IRH depth at the ice core (Δa Δdepth ) and the age uncertainty associated with the ice core at the IRH depth (Δa core ), following where (Δa core ) is a function of the age of the individual IRH at the ice core site (Sigl et al., 2016) and the published uncertainty associated with the ice core age (1% and 3% for ages ranging between 0-15 ka and 15-31 ka BP, respectively; Sigl et al., 2016), while (Δa Δdepth ) is a function of the depth uncertainty of each IRH at the ice-core site. Since the uncertainty in the EM wave through the ice increases with depth, using the maximum uncertainty calculated on the deepest IRH to calculate Δdepth at a catchment scale (see Section 2.3) would result in less accurate age uncertainties at the ice core. We have therefore calculated a depth uncertainty for each individual IRH at the ice core, and undertook the same empirical error analysis to calculate Δdepth at WD2014. This resulted in IRH-specific radar depth uncertainties which we used to calculate the age uncertainty for each IRH at WD2014, as per Equation 2. Whilst Δa comb represents the combined maximum uncertainty from the radar and the ice-core chronology, we found that our IRHs are systematically lower in the ice column compared with strong peaks in acidity concentrations at WD2014 matching closely the age and depth of our IRHs and which we can assume to be the likely cause of our IRHs (see Section 4.2). To account for this offset in ages between the IRHs and the strong sulfate peaks observed at WD2014, we calculated a total age uncertainty (Δa total , Table 3) which represents the maximum age difference between our IRHs and the sulfate peaks at the ice core. This was obtained by adding a systematic factor of 0.22 ka to Δa comb , which represents the total age difference between the maximum IRH age calculated using Δa comb and the age of the strong sulfate peaks (see Section 4.2). We provide the total uncertainty values in Table 3 and Section 3.2.

Age-Depth Modeling
To provide an independent validation of our ice-core derived IRH ages, we also applied the Dansgaard and Johnsen (1969) 1-D vertical ice-strain rate model to derive approximate dating of the IRHs traced over the central PIG catchment. This model has been used previously to date IRHs across West Antarctica (Ashmore et al., 2020;Corr & Vaughan, 2008;Karlsson et al., 2012Karlsson et al., , 2014, assess divide migration (Waddington et al., 2005), and calculate past accumulation rates at or near ice divides (Jacobel & Welch, 2005;Siegert & Payne, 2004). We chose the Dansgaard-Johnsen model here for its simplicity and as it allows us to test the effect of ice deformation on the ages of our IRHs. However, we note that other alternatives exist such as the Nye (1957) and Lliboutry (1979) models, or the more developed quasi-Nye model (MacGregor et al., 2015).
Under the assumption that the ice sheet is, and has been, in steady state, close to an ice divide, the Dansgaard-Johnsen model gives where t (ka; thousand years) is the age of an IRH, H (m) is the ice thickness (assumed constant in time), h (m) is the thickness of the basal shear layer, a (in m a −1 ice-equivalent) is the average accumulation rate since deposition of the IRH, and z (m) is the elevation of the IRH above the bed (Dansgaard & Johnsen, 1969).
For this model, several assumptions are made: (a) negligible horizontal velocity component; (b) time-averaged accumulation rates and no temporal change in accumulation patterns; and (c) constant ice deformation from the surface to some depth, h, below which vertical strain rate is assumed to decrease linearly toward the bed. Considering the above, we initiated the model on the PIG-PASIN data at two sites (A and B in Figure 1) located ∼50 km from the ice divide where horizontal ice flow is minimal (<3 m a −1 ), the ice BODART ET AL.

10.1029/2020JF005927
8 of 20 is thick (>3 km) and the bed relatively flat. Site A (80.15°S, 101.56°W) was selected due to its relative proximity within PIG to WD2014 (∼215 km). At this site, R1-3 were traced, as well as R4. This provided us with initial constraints for age-depth estimates for the upper IRHs (namely R1-3), and allowed us to evaluate the model results based on the approximate known age of R4. To ensure representativeness, however, we also selected a second site, Site B (79.87°S, 100.03°W), where R1-3 were traced but not R4.
We based our estimates for a in the equation on advection-corrected accumulation rates from the WD2014 ice core  for each IRH R1-4, and with current accumulation estimates to correct for any elevation-dependent change in accumulation between the WD2014 site and our PIG Sites A and B. Tentatively treating our R1-3 as broadly equivalent to three of Siegert and Payne's (2004) dated IRHs based on depth associations at three crossovers (see Text S2, Table S3), we derived mean advection-corrected accumulation rates at WD2014 for each reference age: 0.247 ± 0.062 m a −1 (3 ka BP, with BP defined as years before 1950 CE), 0.248 ± 0.062 m a −1 (5 ka BP), and 0.243 ± 0.061 m a −1 (7 ka BP), as well as a rate of 0.226 ± 0.051 m a −1 (17.5 ka BP) based on the intersection with Jacobel and Welch (2005).
The errors correspond to uncertainties in the firn-densification model used by Fudge et al. (2016). These provide us with estimates of what would be required to reproduce each layer if accumulation had remained constant between the time of the deposition of the layer and the present at WD2014. Under the assumption that spatial accumulation patterns have not changed during the Holocene over the WAIS Neumann et al., 2008;Siegert & Payne, 2004), and considering that accumulation rates at the Ice Core are generally smaller than at Sites A and B (Table S4), we use modern accumulation rates from modeled and observational data to calculate the regional difference between accumulation at WD2014 and our  (Neumann et al., 2008) (hereafter referred to as MED14) (Text S2). From these data sets, we calculate a percentage of change between WD2014 and Sites A and B and apply this to the mean advection-corrected rates calculated at WAIS Divide for R1-4 (Table S4). Together, these provided us with a range of realistic values of a for each IRH at Sites A and B to use as input into the 1-D model.
The thickness of the basal shear layer, h, is largely unknown as it is dependent on accurate knowledge of the bed topography and temperature of the ice (Neumann et al., 2008). Previous studies have used a value of  400 m h for Greenland and West Antarctica Jacobel & Welch, 2005;Karlsson et al., 2012;Siegert & Payne, 2004), whilst Karlsson et al. (2014) and Ashmore et al. (2020) explored the effects of fuller ranges of   100 m 1200 m h . We refined this range to 0.2H ≤ h ≥ 0.3H (Neumann et al., 2008), rounding to the nearest 100, hence investigating the effect of h ranging from 700 to 1,100 m at both sites (Text S2). We note, however, that large uncertainties in basal deformation at WD2014 (Cuffey et al., 2016Fudge et al., 2019) could result in h values being smaller than 20% of the ice thickness and thus lead to an overestimation of our ages (see Text S2).

Englacial Stratigraphy
We successfully traced four IRHs R1-4 across a large proportion of the PIG catchment, including in areas where annual velocities reach up to ∼350 m a −1 (Figure 4). The most extensive IRH traced in our study is R2, closely followed by R3 (Figure 4), with mean depths across the catchment of 1,175 and 1,463 m, respectively ( Table 2). The shallowest IRH, R1, was located on average at ∼30% of the ice depth, whilst the deepest, R4, was on average found at ∼68% depth ( Table 2).
The traceability of R1-3 does not vary greatly and is primarily constrained by topography (Figures 4a-4c). By contrast, R4 was only detected across the upper Thwaites/PIG catchments (Figure 4d), even though BODART ET AL.
10.1029/2020JF005927 9 of 20 it has previously been detected much further north into the PIG basin in the ITASE survey (Jacobel & Welch, 2005), likely due to the different frequency range used by the two radar systems. We come back to this point in Section 4.1. We were also able to trace R1-3 in the upper parts of the Institute and Möller ice-stream catchment, and R2-4 in the upper parts of the Thwaites catchment toward the WD2014 site ( Figure 4). The traced IRHs are generally deeper southward of the onset of PIG tributaries 7 and 9 and at the center of the PIG catchment, and relatively shallow at its southern margin and at the divides with Thwaites Glacier and Institute Ice Stream (Figures 4e-4h). We were unable to identify the IRHs in several locations, mainly north of the main trunk of PIG near the Hudson Mountains range and west of tributary 6 ( Figure 4a). We were also unable to detect continuous IRHs in any PIG-PASIN profiles traversing the main trunk and tributaries of Thwaites Glacier, nor those that cover the main trunk and fast-flowing tributaries of PIG (Figure 4).

Age-Depth Estimates
Having clearly identified R2-4 near the WD2014 site, we attempt to date these using the WD2014 chronology. The OIB-MCoRDS2 radar profile passes within ∼1.2 km of the ice-core site, and the stable ice conditions in the area mean that flow-induced disturbance on layer geometry is relatively limited (Laird et al., 2010). Following MacGregor et al. (2015), we calculate the unweighted mean reflection depth within a distance of ±250 m along transect from the trace that is closest to the ice-core site to obtain Δa Δdepth , resulting in mean depths at the ice core of 1,060 ± 7 (R2), 1,430 ± 9 (R3), and 2,371 ± 14 m (R4) ( Table 3). Considering the radar-depth and ice-core uncertainties (Equation 2), and to account for the age offset between our IRHs and the strong sulfate peaks at the ice core (see Sections 2.4.1 and 4.2), we determined the age and associated age uncertainty for each IRH at WD2014 as: 4.72 ± 0.28 (R2), 6.94 ± 0.31 (R3), and 16.50 ± 0.79 ka (R4) ( Table 3).
To compare the ages independently from the WD2014 chronology and provide an approximate age-range estimate for our shallowest isochrone R1, we use the 1-D model at Sites A and B. The age estimates returned from the 1-D model at both sites are as follows: R1 (2.31-2.92), R2 (4.46-5.82), R3 (6.75-9.15), and R4 (19.69-26.87 ka) ( Table 4).
The ages calculated for R2-3 at WD2014 (Table 3) are within the upper and lower bounds of the modeled age-range estimates from the 1-D model ( Note. We provide these for both depths below the surface and depth as a fraction of ice thickness. "1σ" refers to one standard deviation, "Range" refers to the minimum and maximum values, and "IQ range" refers to the interquartile range (25th and 75th percentile). A maximum uncertainty of ±17 m is assumed here. Note. Column "a (ka)" refers to the IRH age obtained from the radardepth and the depth at the WD2014 ice core. Column "Δa comb " refers to the combined age uncertainty from the radar and the ice-core chronology, whilst "Δa total " refers to the maximum age uncertainty of our IRHs calculated from the age difference between our IRHs and the strong sulfate peaks at WD2014 (see Sections 2.4.1 and 4.2). Abbreviation: IRH, internal reflecting horizon. 4. Discussion Karlsson et al. (2014) traced two distinctive IRHs through the middle ice depths across parts of the central PIG catchment using the same PIG-PASIN data set as that used here, but only focusing on flight lines flown at constant elevation and only exploiting the data in its chirp mode. This earlier study highlighted the existence of a distinctive IRH package between an upper bound, "Layer 1", approximately dated to 5.3-6.2 ka, and a lower bound "Layer 2", approximately dated to 8.6-13.4 ka. Here, by additionally exploiting the full spatial extent of the PIG-PASIN data set, the simultaneously acquired pulse-mode PASIN data, and complementing these with recent OIB-MCoRDS2 data, we have expanded the reach of that earlier radiostratigraphy across the fuller PIG catchment, and across the ice divides into neighboring regions, notably Thwaites Glacier and Institute Ice Stream. Direct comparison between both sets of results suggests that Karlsson et al.'s (2014) Layers 1 and 2 are equivalent to the IRHs traced in this study as R2 and R3, with a median difference ranging between 6 and 12 m, which is within the depth uncertainty of the IRHs ( Figure S3).

Table 4 Modeled IRH Age-Range Estimates (ka) Returned From the 1-D Steady-State Model for Varying Accumulation Data Sets (see Section 2.4.2) and Basal Shear Layer Thickness (h, in m) Scenarios at Sites A and B for IRHs R1-4 (see Section 2.4.2)
et al. 2020), and hence provides additional evidence that we observe the same IRHs across both catchments. Two sets of parallel profiles, laterally offset by ∼1.5 km, and acquired across the PIG/Institute Ice Stream divide in the PIG-PASIN and IMAFI-PASIN data sets (Figure 1), provide a further opportunity to confirm these equivalences with data from the same radar system. Only in three short sections of these transects could we compare our IRHs with those from the IMAFI-PASIN study (inset Figure S3a); in these locations, BODART ET AL.

10.1029/2020JF005927
12 of 20 Also shown are the IRH ages (ka) (see Section 3.2) for R1 (age-range estimate from 1-D model) and R2-4 (ages from WD2014 ice-core intersection). For (a-d), lower (blue) values correspond to relatively deep IRH depths, higher (yellow) values correspond to shallow IRH depths. Background is bed elevation in meters (referenced to the WGS84 ellipsoid) from BedMachine (Morlighem et al., 2020). For (e-h), lower (yellow) values correspond to the shallowest IRHs, higher (purple) values correspond to the deepest IRHs. Background is ice thickness in meters from BedMachine (Morlighem et al., 2020). we could not identify R1 and R3. Nevertheless, at two intersections (black arrows in inset on Figure S3a), the respective depths for PIG-PASIN R2 and IMAFI-PASIN H2 were 794 and 797 m at Intersection 1 and 776 and 778 m at Intersection 2, respectively, which is remarkably close considering ice thickness in this area exceeds 2 km. This, alongside the crossovers on the OIB-MCoRDS2 data, gives us high confidence that our R2, Ashmore et al.'s (2020) H2, and therefore Karlsson et al.'s (2014) Layer 1, all represent the same internal marker in the ice. This study, by using additional data that allowed direct dating at the WD2014 site, is therefore able to ascribe more accurate and precise ages to the IRH package ranging across PIG and Institute and Möller Ice Streams of 4.72 ± 0.28 ka (Layer 1/H2/R2) and 6.94 ± 0.31 ka (Layer 2/H3/R3), respectively, based on the WD2014 ice-core chronology.
We also note that all three studies identify R2 as their most spatially extensive IRH, indicating the presence of a particularly ubiquitous isochrone, similar in age to a 4.72 ± 0.24 ka isochrone detected and also extensively mapped elsewhere across central West Antarctica ). Whilst we were not able to provide a more refined age to our shallowest IRH, R1, from direct intersection of the WD2014 Ice Core, the 1-D model returned an age-range estimate (2.31-2.92 ka) that is in broad agreement with that of Ashmore et al. (2020) (1.9-3.2 ka; their H1) and Siegert and Payne (2004) (3.10 ± 0.16 ka; their L07). Together, these studies demonstrate considerable promise for unifying an age-depth stratigraphy across the WAIS back to at least ∼7 ka, while tying our IRHs to the WD2014 Ice Core has yielded more accurate, and younger, ages, for the isochrones detected across PIG and, by extension, Institute and Möller Ice Streams.
The age assigned to R4 at WD2014 (16.50 ± 0.79 ka) is slightly younger than the 17.5 ± 0.5 ka layer tied by Jacobel and Welch (2005) to the Byrd Ice Core (Hammer et al., 1997), although there is an overlap of 0.29 ka when fully accounting for the age uncertainties. We offer two potential explanations for this disparity. First, the low-frequency ground-radar system used as part of the ITASE survey has a much longer wavelength than the high-frequency airborne systems used here, meaning that the 17.5 ± 0.5 ka layer appears as a single-amplitude peak measuring tens of meters in thickness (cf. Jacobel & Welch, 2005), whereas the shorter-wavelength on the airborne radars allows for the delineation of individual peaks, thus resolving the strong singular reflector from Jacobel and Welch (2005) as a series of closely spaced reflectors. As a result, when attempting to connect the ITASE profile with the airborne radar data, it is likely that the closest bright reflector identified on the airborne radar forms the upper part of the wider reflector imaged by Jacobel and Welch (2005), thus leading to younger ages at the intersection with the WD2014 Ice Core. Second, the uncertainties in the radar data at the intersection between OIB-MCoRDS2 (±14 m) and Jacobel and Welch's (2005) profile (±10 m) increase the chance to misinterpret the correct position of the 17.5 ka layer over the airborne data, although we show in Table S2 that the mean depth difference between R4 and Jacobel and Welch's (2005) layer is < 18 m, which is within the uncertainty range of both studies. Whilst these points are relevant when comparing the ages of R4 at WD2014 with the age of Jacobel and Welch's (2005) layer, it is worth mentioning that the exact age and depth of the strong reflector at WD2014 are known from electrical conductivity and chemistry measurements. At the ice core, this layer is characterized by nine distinctive peaks ranging in depths between 2,420 and 2,427 m and dated at 17.75 ± 0.19 ka (McConnell et al., 2017;Sigl et al., 2016), a full 35 m below the depth of R4 at WD2014. Even taking into account the maximum depth of our IRH along the ±250 m transect (2,378 ± 14 m; see Section 3.1), R4 is still found 28 m above the depth of the 17.75 ± 0.19 ka at WD2014. Considering all the above, it is likely that R4 is not the same layer as the strong volcanic layer dated at 17.75 ± 0.19 ka at WD2014 (McConnell et al., 2017), but rather forms the upper part of the wide reflector imaged by Jacobel and Welch (2005) in the ground-radar data.

Linkage with the WAIS Divide Ice-Core Record
Whilst determining the cause of R4 remains ambiguous due to the limitations mentioned above, the existence of R2 and R3 offers an opportunity to link them directly to the ice-core sulfate record at WD2014. High sulfate content from volcanic sulfuric acid is known to correspond to high acidity levels in englacial layers in ice cores (Castellano et al., 2005;Gow & Williamson, 1971;Hammer et al., 1997;Millar, 1982) and, because the radar is sensitive to acidity contrasts (Fujita et al., 1999;Millar, 1981), we can attempt to link the sulfate record at the ice core with our IRH stratigraphy. Figure S4 shows the presence of three large peaks in sulfate concentration at the WD2014 ice core, which are particularly close in age and depth to IRHs R2-3 traced on the OIB-MCoRDS2 profile near WD2014. In particular, a layer dated at 4.94 ka (depth: 1,099 m) BODART ET AL.

10.1029/2020JF005927
13 of 20 contains sulfate concentrations that are unmatched (405 μg/kg) for much of the core up until a depth of ∼2,400 m (equal to the last ∼18,000 years BP) ( Figure S4). Even taking into account the entire profile, this layer contains the fourth largest amount of sulfate concentrations in the last ∼68,000 years BP. We also notice the presence of two closely spaced peaks in the sulfate record which are dated at 7.25 ka (depth: 1,475 m; sulfate concentration: 306 μg/kg) and 7.64 ka (depth: 1,526 m; sulfate concentration: 271 μg/ kg), corresponding to the 9th and 10th highest sulfate concentrations on record ( Figure S4b). Not only do these ages match closely the age of R3 at the ice core, they also match the characteristics of R3, which is often found as a couplet across most of Pine Island, upper Thwaites, and Institute and Möller ice-stream catchments on the airborne radar data (Figures 2 and S1). Additionally, the second largest peak on record before ∼18,000 years BP is found at a depth of 584 m and dated at 2.45 ka (sulfate concentration: 309 μg/ kg), which falls within the modeled age-range estimate for R1 (2.31-2.92 ka) at Sites A and B (Table 4, Figure S4a).
Whilst this offers us the opportunity to directly link our IRHs to the WAIS Divide record, we note that the depths of R2-3 at WD2014 are slightly shallower (R2: 1,060 ± 7 m; R3: 1,430 ± 9 m) than the sulfate peaks in Figure S4, resulting in slightly younger ages at the ice core. We cannot exclude the possibility that we traced a layer that is slightly above R2 and R3 at the Ice Core, although this is unlikely as we base our tracing on depth intersections ( Figure S2) and IRH characteristics ( Figure S1). Even taking into account the maximum depth of R2-3 along our ±250 m transect to account for the fact the OIB-MCoRDS2 line did not fly directly over the WD2014 site but instead ∼1.2 km away (see Section 2.4.1), R2 (1,069 ± 7 m) and R3 (1,438 ± 9 m) would still be found 23 and 28 m higher than the sulfate peaks at the ice core, respectively. Whilst this is a relatively small disparity considering ice thickness in the area exceeds ∼3.5 km and that we are effectively comparing airborne-radar data (meter-scale accuracy) with ice-core data (mm-scale accuracy), the reason for our IRHs not aligning more closely with the sulfate peaks remains unclear. One potential explanation could relate to the distance between our transect and the location of the WD2014 ice-core site. Although Laird et al. (2010) suggested that flow-induced disturbance on layer geometry is limited in the area around the WD2014 site, changes in bed roughness were found to affect englacial stratigraphy near WD2014. This could lead to small undulation in IRH elevations between our transect and WD2014 and thus cause in several meters of discrepancy. To acknowledge this, and considering that the sulfate peaks are most likely the cause of our IRHs as we show above, we have increased the age uncertainty of our IRHs to account for the offset between our IRH ages and the age of the sulfate peaks (see Section 2.4.1, Table 3). This results in more conservative uncertainties for our deeper three IRHs dated at the ice core: 4.72 ± 0.28 (R2), 6.94 ± 0.31 (R3), and 16.50 ± 0.79 ka (R4).
By linking three of our four IRHs to the sulfate record at WAIS Divide, we can hypothesize that the origin of our spatially extensive IRHs is from past explosive volcanic activity during the Holocene. Previous studies in Antarctica have demonstrated the correspondence between bright reflectors in radar data and past volcanic activity (e.g., Corr & Vaughan, 2008;Jacobel & Welch, 2005). Karlsson et al. (2014) previously attempted to link their deeper layer (Layer 2/R3) to acidity peaks at Byrd Ice Core; however, the absence of a direct link between the PIG catchment and a complete ice-core chronology was lacking at the time. The evidence presented here suggests that our IRHs may also originate from past explosive volcanism; but, the precise source of these eruptions, whether regional or global, remains unknown.

Accumulation Rate and IRH-Age Comparison
The correspondence in isochrone-age estimates for IRHs R2-3 derived from intersecting the WD2014 site (Table 3) and using the 1-D model (Table 4) at the PIG/Thwaites divide (∼250 km away) (our Sites A and B; Figure 1) suggests that accumulation patterns have remained broadly similar across the Amundsen-Ross divide for at least the last ∼7 ka. Whilst this is based on a relatively limited amount of data points, it complements previous studies Koutnik et al., 2016;Neumann et al., 2008), including Siegert and Payne (2004) who, using the same SPRI/NSF/TUD radar transect as that in Figure 1, concluded that accumulation patterns have remained stable over the last 6.4 ka. We suggest future research make use of the accurately dated IRHs provided here to model Holocene accumulation rates and patterns, as well as regional ice-sheet balance velocities, as previously conducted over Greenland (e.g., MacGregor et al., 2016) BODART ET AL.

10.1029/2020JF005927
14 of 20 and on individual sections of the WAIS Neuman et al., 2008). This will provide additional information on the terrestrial ice-sheet history of the ASE during the Holocene, and in turn help us to constrain better the future of the WAIS.
Previous studies have successfully combined ice-core records with modeled modern-day accumulation rates to reconstruct Holocene accumulation (Cavitte et al., 2018;Fudge et al., 2016;Nielsen et al., 2018), although non-climatic noise in the observations and model biases have resulted in small discrepancies between icecore and model reconstructions (Cavitte et al., 2020;Dalaiden et al., 2020). When assessing the ability of the 1-D model to reproduce the ages for R2-3 derived at the WD2014 Ice Core, we find that the best match (to within < 10%) is achieved using the modern accumulation rates provided by the MED14 and RACMO2 products. This is not surprising as both have higher spatial resolution than MAR and ART06, but it also likely reflects the fact that MED14 is an observational product and that RACMO2 has been shown to agree well with geophysical estimates of accumulation rates (Lenaerts et al., 2012;Medley et al., 2014;van Wessem et al., 2018;Wang et al., 2016). In contrast, when using present-day accumulation estimates from ART06 and MAR to calculate past accumulation rates, model-derived ages are up to 1.1 ka (∼23%) greater for R2 and 2.2 ka (∼32%) greater for R3 compared with ice-core derived ages (Tables 3 and 4). This discrepancy is primarily dominated by different modern accumulation gradients estimated between WD2014 and the PIG/ Thwaites divide (i.e., Sites A and B), with the MED14 and RACMO2 products suggesting a slightly more homogenous gradient than ART06 and MAR (Table S4). Lower in the ice, the poor correspondence between the age of R4 derived by links to the WD2014 (16.50 ± 0.79 ka) relative to the age returned by the 1-D model (19.69-26.87 ka) is worthy of investigation. Even taking into account the maximum age uncertainty at the ice core, the minimum and maximum age returned by the 1-D model is 2.6 (15%) and 9.8 ka (57%) greater than at the ice core (Tables 3 and 4), a difference that cannot solely be attributed to the different modern-day accumulation gradients mentioned above. The most likely explanation is that the assumptions required for the 1-D model (see Section 2.4.2) break down for older IRHs, where local accumulation rate is no longer a primary factor in determining the depth of an IRH. This could be due to complex flow dynamics such as longitudinal strain or lateral shearing at the boundary between slow and fast-flowing ice, resulting in high internal stress impacting IRH stratigraphy in the deeper part of the ice column (Waddington et al., 2007). Moreover, R4 (16.50 ± 0.79 ka) was deposited pre-Holocene as the WAIS was transitioning from a glacial to an interglacial period during which ice thickness has likely not remained constant (Golledge et al., 2014;Johnson et al., 2017), implying possible changes in ice-flow configurations for which the steady-state model is not able to account.

Characteristics of Englacial Stratigraphy
Previous research over East Antarctica has shown that common bright reflectors can be interchangeably traced over long distances using radar systems operating at different center frequencies (Cavitte et al., 2016;Winter et al., 2017). Our findings provide further evidence of this over West Antarctica, having successfully identified common IRHs across different airborne radar systems. However, although IRHs younger than 7 ka can be traced widely across the WAIS using existing data sets, tracing deeper, pre-Holocene IRHs has not been widely possible across PIG (this study) nor the Weddell Sea Sector (Ashmore et al., 2020). Relative to the interior of East Antarctica, where much lower snow accumulation and ice-flow velocities have facilitated the tracing of isochrones pre-dating the LGM (∼20 ka BP) and even the past glacial-interglacial periods (up to ∼366 ka BP) (Cavitte et al., 2016;Parrenin et al., 2017;Steinhage et al., 2013;Winter et al., 2019), the extremely variable deep-ice conditions in the WAIS will challenge the recovery of pre-Holocene radiostratigraphy. Compounding the challenge, Ross et al. (2020) have demonstrated that large packages of ice older than ∼16 ka in the Weddell Sea sector of the WAIS are rheologically different to the ice above, containing large proportions of deformed and folded ice. These packages typically show poor continuity of englacial stratigraphy across Institute and Möller Ice Streams (Bingham et al., 2015) and, indeed, where we could see IRHs deeper than R4 in PASIN and MCoRDS2 for this study, very few were continuous for long distances. Over other parts of the WAIS, an IRH dating back to 24.9 ± 0.3 ka has been traced in limited radar profiles connecting the Byrd and WAIS divide ice cores, where it was found at 68% and 80% of ice depth at Byrd and WD2014, respectively ; however, they were also unable to recover deeper continuous IRHs more widely.

10.1029/2020JF005927
15 of 20 Overall, with the existing data sets available across the WAIS, the prospects for tracing and dating Holocene radiostratigraphy widely across the ice sheet with existing data are excellent, but diminish rapidly for older ice, going back to the LGM and beyond. Yet, much deeper, and thus older IRHs, are visible throughout the ice column with ground-based radars (e.g., Bingham et al., 2017;King, 2011;Laird et al., 2010) and hence the interrogation of older ice in the WAIS may be best suited to strategic ground campaigns that can be linked into the airborne-derived radiostratigraphy. In the PIG catchment, older ice is suggested by our results to lie below the PIG/Thwaites divide, where on average ∼900 m of ice (30% of the mean ice thickness) underlies R4 (∼17 ka) ( Figure S5).

Conclusion
We have identified four spatially extensive IRHs in airborne radar surveys that are present across much of the Pine Island Glacier catchment in West Antarctica. Extending into neighboring Thwaites Glacier and Institute Ice Stream, these IRHs can be considered isochrones that span the late Pleistocene and Holocene, with ages of 2.31-2.92, 4.72 ± 0.28, 6.94 ± 0.31, and 16.50 ± 0.79 ka derived from intersecting the WAIS Divide Ice Core and the use of a 1-D ice-flow model. Our most spatially extensive IRH, R2, is remarkably similar in age and depth to another extensive IRH previously identified by other studies over Pine Island Glacier, Institute and Möller Ice Streams, and the Marie Byrd Land region. More broadly, we have also shown that our IRH package is similar to previously traced IRHs over the Weddell Sea sector of the WAIS, which, together with the Pine Island Glacier catchment represents ∼20% of West Antarctica. Finally, we have shown that our upper three IRHs correspond to large peaks in sulfate concentrations at the WAIS Divide Ice Core, suggesting that our IRHs are of volcanic origin.
When assessing the presence of older ice across the catchment, we observe that the relative proportion of ice older than R4 in the ice column is limited and does not contain many continuous reflections. Indeed, we find that the deepest (and thus oldest) continuous IRH identified in this study, R4, is found at an average depth of 68% in the ice column despite its age (∼17 ka), only representing 25% of the estimated age of the oldest ice recovered at the WAIS Divide Ice Core (∼68 ka). This indicates that the majority of ice older than the LGM is found within the bottom ∼30% of the ice thickness across PIG/Upper Thwaites. Whilst this is to be expected as the age-depth profile of an ice sheet does not increase linearly, the absence of continuous reflections dating back to the LGM and older currently limits our ability to reconstruct longer-term changes using existing airborne data sets.
As isochronous features, the dated IRHs generated here offer a new set of large-scale boundary conditions that could be a valuable resource, if incorporated into ice-flow models seeking to improve our understanding of past ice-sheet evolution. We anticipate that these well-dated IRHs will provide constraints for models simulating past accumulation rates and patterns, which in turn will shed more light onto the terrestrial ice sheet history of this very sensitive catchment of the WAIS.