A 3 m.y. Record of Volcanism and Glaciation in Northern British Columbia, Canada

The Tuya-Kawdy region of northern British Columbia is well established as a place where glaciation and volcanism overlapped in space. However, no modern work has integrated observations from the region’s volcanic and glacial deposits with geochronologic constraints to summarize how they might overlap in time. Here, we provide a general overview of such characteristics and 23 new 40Ar/39Ar eruption ages of glaciovolcanic deposits ranging from 4.3 Ma to 63 ka to constrain the timing, location, and minimum thicknesses and distributions of coincident ice. Subaerial lava fields interspersed with glaciovolcanism record periods of ice-sheet absence in presumably warmer climate conditions. These generally coincide with interglacial marine isotope stages. Many of the volcanoes have a secondary record of posteruption glacial modification, cirques, erratics, and mega-lineations, which document later climate changes up to the present. We used edifice-based terrain analysis to reconstruct changes to local minimum Cordilleran ice-sheet thicknesses, extents, and flow directions at specific locations and times during the late Pliocene and the Pleistocene. *published posthumously Edwards, B.R., Russell, J.K., Jicha, B., Singer, B.S., Dunnington, G., and Jansen, R., 2020, A 3 m.y. record of volcanism and glaciation in northern British Columbia, Canada, in Waitt, R.B., Thackray, G.D., and Gillespie, A.R., eds., Untangling the Quaternary Period—A Legacy of Stephen C. Porter: Geological Society of America Special Paper 548, p. 231–257, https://doi.org/10.1130/2020.2548(12). © 2020 The Authors. Gold Open Access: This chapter is published under the terms of the CC-BY license and is available open access on www.gsapubs.org. OPEN ACCESS GO LD Downloaded from http://pubs.geoscienceworld.org/books/book/chapter-pdf/5325010/spe548-12.pdf by Dickinson College user on 03 August 2021


INTRODUCTION
Glaciovolcanism occurs wherever volcanoes and glaciers coincide in space and time (Smellie, 2000;Kelman et al., 2002;Russell et al., 2014;Edwards et al., 2015;Smellie and Edwards, 2016). Their interactions leave behind distinctive deposits like pillow lavas and palagonitized tephra, uncommon lithofacies associations, and unique landforms such as tuyas that record paleoenvironmental conditions during eruptions. The products of such glaciovolcanism may be used to map and date ice distributions and the waxing and waning of glaciations.
The coincidence of volcanism and glaciation has long been recognized in western Canada ( Fig. 1; Kerr, 1926;Mathews, 1947Mathews, , 1951Mathews, , 1952Grove, 1974;Souther et al., 1984;Edwards and Russell, 1994;Spooner et al., 1995;Hickson et al., 1995;Hickson, 2000;Edwards et al., , 2009Kelman et al., 2002;Wilson and Russell, 2018). Previous workers have inferred that the Tuya-Kawdy area of northwestern British Columbia, Canada, in particular has been inundated frequently by major ice sheets over the past ~3 m.y. (Ryder and Maynard, 1991;Barendregt and Irving, 1998;Clague and Ward, 2011;Edwards et al., 2011). During this same time, tectonic forces have caused mafic alkaline volcanism within the overlapping northern Cordilleran volcanic province (Edwards and Russell, 2000) to produce a landscape with at least 35 glaciovolcanoes and a few subaerial volcanoes. Here, we used this volcanic record to reconstruct the glacial history of the Tuya-Kawdy region of British Columbia. We examined historic regional mapping studies (Kerr, 1926;Watson and Mathews, 1944;Gabrielse, 1970Gabrielse, , 1998, documented volcanic lithofacies through our own field work, analyzed the morphology of individual edifices, and integrated glacial geomorphology and deposits with new volcanic geochronology.   (Tuya-Kawdy area), and the dotted outline shows the map extent for Figure 14. The extents of the Cordilleran and Laurentide ice sheets, the ice-free zones, and the inferred directions of ice movement are from Clague (1989). Star-Vancouver, British Columbia, Canada. Inset map shows the broader context of northwestern North America including the approximate modern-day continental drainage divides and three major rivers (Stikine, Mackenzie, Yukon) that flow into three different marine environments (Pacific and Arctic Oceans, Bering Sea). Each of the three rivers gets water from within the study area.

Paci c Plate
The Tuya-Kawdy area of northern British Columbia is a classic location for tuya research (Fig. 1); for example, see the seminal work of Mathews (1947). "Tuya" is a word in the native Tahl'tan language meaning "baby water," and it is an element of many local geographic names (e.g., Tuya Lake, Tuya Butte, and Tuya River). Tuya Butte, which is a large, distinctive flat-topped mountain ( Fig. 2; Mathews, 1947), is the namesake for all other volcanoes characterized by morphologies that are a direct result of eruptions under ice (elaborated below).

REGIONAL SETTING OF VOLCANISM IN WESTERN CANADA
Western Canada was characterized by sporadic volcanism for much of the Cenozoic (Souther and Yorath, 1991). Until ca. 20 Ma, the North American-Kula/Farallon subduction zone caused volcanism along most of the coast of British Columbia and into Yukon Territory (Edwards and Russell, 1999;Madsen et al., 2006). The Garibaldi belt, the northernmost extension of  Table 1. Dashed elevation contours for 1000 and 1500 m above sea level are included for reference; the 1000 m contours are labeled, and the 1500 m contours not. The dashed-white box shows approximate area of Figure 10A.
the Cascade volcanic arc into southwestern British Columbia, is the only remnant of that extensive subduction still active today (Fig. 1). After a hiatus from subduction magmatism, extensionrelated alkaline volcanism began around 10 Ma in several parts of British Columbia (Souther and Yorath, 1991). The alkaline magmatism has been attributed to (1) the formation of a slab window beneath northern British Columbia (Thorkelson and Taylor, 1989;Madsen et al., 2006), (2) transtension along the Queen Charlotte strike-slip boundary that separates the Pacific and North American plates (Souther and Yorath, 1991;Souther, 1992;Russell, 1999, 2000), and (3) extension due to lithospheric relaxation after an extended period of terrane accretion and subduction (Manthei et al., 2010). Volcanism since 3 Ma may have been affected by isostatic adjustments because of glacial loading and unloading (Grove, 1974;Edwards and Russell, 1999).
The volcanism has been divided into a number of regions that span British Columbia and Yukon Territory southeast to northwest ( Fig. 1), including the Wells Grey-Clearwater volcanic field Hickson and Vigouroux, 2014), the Anahim volcanic belt (Bevier et al., 1979;Kuehn et al., 2015), the Garibaldi belt of the Cascade volcanic arc (Green, 1990;Kelman et al., 2002;Wilson and Russell, 2018), and the northern Cordilleran volcanic province (Edwards and Russell, 2000). All have glaciovolcanic deposits recording volcanism coincident with climatic switches that caused localized glaciers and the regional Cordilleran ice sheet. The history of the Cordilleran ice sheet has been partly tracked (Ryder and Maynard, 1991;Barendregt and Irving, 1998;Froese et al., 2000;Clague and Ward, 2011).

TUYA-KAWDY VOLCANIC FIELD
The northern Cordilleran volcanic province hosts the Tuya-Kawdy plateau, where W.H. Mathews (1947) first recognized a distinctive volcanic landform, which he named "tuya," formed from sustained eruptions from a single vent beneath a continental ice sheet. The Tuya-Kawdy volcanic field is located within a region consisting of large plateaus and small mountain ranges that together make up parts of the Kawdy Plateau and the Cassiar Mountains (Figs. 2 and 3). Most peaks and ridges within the mountains are separated from each other by large U-shaped valleys.
The physiography of the region varies significantly. To the southwest, the Kawdy Plateau rises above surrounding river valleys to a basal elevation of ~1300 m above sea level (a.s.l.), upon which six volcanic mountains have been constructed (Fig. 2). The volcanoes are the highest points on the plateau (Fig. 3A). To the north and northwest, various isolated ranges together form the Cassiar Mountains, where peaks reach above 2100 m. Many of the volcanoes are located within the valleys and between the higher peaks of the Cassiar Mountains (Fig. 3B). Kerr (1926) and Watson and Mathews (1944) initially mapped parts of this region, and Gabrielse (1970Gabrielse ( , 1998 finished the modern regional mapping as part of a 1:250,000 project. Through these works, the "Tuya Formation" was defined as a regionally extensive stratigraphic unit (Kerr, 1926). Kerr (1926) recognized that these volcanic deposits, which he referred to as part of the "Tuya lava field," probably erupted within the Pleistocene, because they had experienced variable amounts of glacial modification. He also speculated that in at least one area, Canyon Creek, a subaerial lava erupted when glaciers were not present. The region-scale mapping greatly facilitated later topical studies (e.g., Allen et al., 1982;Moore et al., 1995;Dixon et al., 2002;Wetherell et al., 2005;Simpson et al., 2006;Edwards et al., 2011;Russell et al., 2013;Ryane et al., 2011;Turnbull et al., 2016). Mathews (1947) not only connected these unique volcanoes to past regional ice sheets but also described a basic model that is still envisioned today for their formation within intraglacial lakes. Mathews enlarged upon Nielsen's (1936) discussion of an eruption from Grimsvötn volcano in Iceland, where initial explosive eruptions occurred within the lake and then growth of an edifice that sometimes breeched the lake surface. Mathews (1947) pointed out that the volcanoes built large deltaic structures into these intraglacial lakes, and that the transition from subaqueous to subaerial lava flows is a stratigraphic marker recording the former lake level. He speculated that tuya structures were more likely to form via eruptions into "intraglacial" lakes (totally confined by ice) than in periglacial lakes. He also noted that the subglacial topography might have played a role in determining the height of lake surfaces and when and how the eruption-formed intraglacial lakes might have drained through "rocky spillways." While a general relation between glaciation and volcanism had been made in Iceland a few decades earlier (Peacock, 1926), the 1947 paper by Mathews was the first widely distributed, detailed explanation of tuya formation and is still worth reading today.
The most recent regional geologic maps from Gabrielse (1970Gabrielse ( , 1998 show >50 deposits that are part of the stratigraphic unit designated as the basaltic Tuya Formation; we confirmed that at least 35 of these are eruptive centers (Table 1). All 35 are known to be basaltic in composition, though detailed, published petrological studies are few (see Table 1). Little is known about most of these volcanoes beyond their general assignment to the Tuya Formation. Results presented here are based on our field studies between 1995 and 2014 documenting the extents and character of the Tuya Formation (Wetherell et al., 2005;Simpson et al., 2006;Edwards et al., 2011;Ryane et al., 2013;Russell et al., 2013;Turnbull et al., 2016). We have visited all 35 distinct edifices, which are dominantly glaciovolcanic in origin.

Tuya Morphologies
Tuya morphology has largely been used to better constrain eruption environments (Jakobsson and Gudmundsson, 2008;Smellie, 2009Smellie, , 2013Russell et al., 2014;Pedersen and Grosse, 2014). Here, we used the classification of Russell et al. (2014) to catalogue variations in tuyas morphologies within the Tuya-Kawdy volcanic field. This classification directly derives from close reading of Mathews (1947), where the word "tuya" was    first applied to volcanoes thought to be directly shaped by eruption into ice. Mathews (1947) emphasized the characteristic flat top of Tuya Butte as one of the defining features of tuyas, but he referred to other volcanoes in the Tuya-Kawdy area (e.g., Kima 'Kho) with more complex shapes as "tuyas" as well. The nomenclature used by Russell et al. (2014), which applies only to volcanoes and not to regionally extensive deposits (e.g., moberg and pillow sheets of Iceland [Jakobsson and Gudmundsson, 2008]), divides tuyas into four types by primary morphology: (1) tindars or linear tuyas (both terms acceptable in Russell et al. [2014], where "tindar" is a long-established morphological term [Jones, 1968]), (2) conical tuyas (called "subglacial mounds" by Hickson, 2000), (3) flat-topped tuyas (or flat-top tuyas), and (4) compound tuyas. Tindars, which are elongate volcanic ridges, were first described by Jones (1968Jones ( , 1970 in Iceland and form during subglacial fissure eruptions. Conical tuyas typically have shapes that are reminiscent of subaerial volcanic cones, except they are extensively palagonitized. Flat-top tuyas have the "classic" tuya morphology, with a flat top formed by subhorizontal, subaerial lavas overlying subaqueous deposits. Compound tuyas consist of single edifices with combinations of morphologies 1-3 and have only been described from the Tuya-Kawdy area. All four types are present in the study area ( Fig. 4; Table 2): Six are flat-top, eleven are conical, two are linear, and four are compound. The ability to recognize the primary (e.g., not altered by later glaciation) morphology of a tuya is useful because it is direct evidence of the broad-scale, glacio-hydrological conditions during the eruption (Smellie and Skilling, 1994;Smellie, 2006;Russell et al., 2014;Smellie and Edwards, 2016). The Russell et al. (2014) classification simplifies an already complex terminology because all glaciovolcanoes, regardless of morphology, are included in the term "tuya."

Diagnostic Glaciovolcanic Lithofacies
Discrimination between glaciovolcanic and subaerial volcanoes may also be achieved by recording the combinations of unique coherent (e.g., lava) and fragmental lithofacies produced in ice-dominated environments (cf. Smellie and Edwards, 2016). While most of these have been found in the Tuya-Kawdy area ( Table 2; Figs. 5 and 6), the most common coherent lithofacies include pillow lavas (Figs. 5A-5D), subaerial lava sheets, and lava masses featuring pervasive, intense, and irregular to radial jointing patterns and joint orientations (Figs. 5E and 5F). Pillow lavas clearly indicate eruptions beneath or into water (e.g., Walker, 1992), including marine settings, lakes, and rivers. Thus, pillow lava alone does not indicate glaciovolcanism. The physiographic, topographic, and stratigraphic context of the pillow lavas all together diagnose glaciovolcanism. For example, Mathews (1947) noted that Ash Mountain has pillow lavas at its base but also within the capping tephra deposits near its summit, 800 m above its base. At 16 glaciovolcanoes in our study area, pillow lava sequences occur where no present-day means of water impoundment exist.
Certain patterns of columnar jointing in lava masses also provide evidence of accelerated and variable cooling history and commonly reveal the presence of water or ice (Long and Wood, 1986;Spörli and Rowland, 2006;Forbes et al., 2014). Two particular jointing characteristics are frequently cited as diagnostic: (1) joints that have irregular, subvertical orientations but that are not parts of colonnade-entablature sets, and (2) joints with unusually small (less than ~1 m) cross-sectional diameters (Lodge and Lescinsky, 2009). At Mathews tuya, massive lavas on its north flank have tube-like morphologies with fractures that converge from the outside to a central plane, forming radial patterns in cross section ( Fig. 5E; Edwards et al., 2011). In other locations, the presence of pervasive subhorizontal jointing indicates lava-water-ice interactions ( Fig. 5F). At least two centers (Klinkit Creek and Mount Josephine) have massive lava emplaced directly on top of pillow lava (Fig. 5G). These may have formed from emplacement of lava into ice tunnels on the volcano flanks that were periodically filled with meltwater, but that later filled with massive lava after the meltwater drained (e.g., Oddsson et al., 2016).
Fragmental lithofacies diagnostic of glaciovolcanism are found at all of the volcanic centers in the Tuya-Kawdy volcanic field (Table 2; Fig. 6). The fragmental lithofacies themselves include conventional pyroclastic deposits and processes (e.g., tuff, lapilli tuff, tuff-breccia; Table 2), but several support a glaciovolcanic origin. Locally abundant tuff and lapilli tuff vary from unstructured to showing moderate size sorting and bedding (Figs. 6A and 6B). In thin section, these tuffs commonly have subequal amounts of blocky and shard-like shapes, which for basaltic eruptions suggests involvement of (melt-)water. Cored armored lapilli also suggest phreatomagmatic activity (Fig. 6C). So do bomb sags (Fig. 6D), which require a "wet" eruption plume or water-saturated tephra. Coarser-grained volcaniclastic deposits commonly contain fragments of pillow lava (Fig. 6E) or large rounded exotic lithic clasts. The pillow fragments derive from secondary explosions disrupting and dispersing pillow lava, and the exotic lithics are potentially derived from melting out of enclosing ice (Fig. 6F). Lastly, most fragmental deposits have a distinctive yellow-orange color ( Fig. 6), which invariably indicates partial to pervasive palagonitization. Palagonite forms by hydration and alteration of vitric basalt (sideromelane) in warm, water-rich environments (e.g., Jakobsson, 1978;Jercinovic et al., 1990;Stroncik and Schmincke, 2002). Though palagonite is not uniquely diagnostic of glaciovolcanism, glaciovolcanic edifices almost always contain partially to extensively palagonitized deposits (Harder and Russell, 2007;Smellie and Edwards, 2016). Its pervasive presence from base to the summit of volcanoes that sit hundreds of meters above the surrounding landscape indicates that the deposits formed in ice-confined, water-dominated environments.

Tuya Ages
Before this study, only one volcano from the area, Mathews tuya, had a published eruption age . Here,    Table 1). Lp-pillow lava, Lm-massive lava, Ld-dike, VB-volcanic breccia, LT-lapilli tuff, T-tuff, D-diamictite, Dg-glacial diamictite. we report 23 new 40 Ar/ 39 Ar ages for volcanoes within the Tuya-Kawdy volcanic field (Table 3; Fig. 7). The ages are calculated relative to 28.201 Ma for the Fish Canyon sanidine (Kuiper et al., 2008) using decay constants of Min et al. (2000). The 40 Ar/ 39 Ar methods are from Jicha et al. (2012). Of these 23 new dates, 18 are tuyas, and five are from deposits interpreted to have erupted subaerially. The ages range from 4.30 Ma to 0.063 Ma (Table 3). The oldest deposits are on the Kawdy Plateau in the southwestern part of the study area, and the other ages are interspersed to the east and north. Of the 18 newly dated glaciovolcanic deposits, eight are flat-topped or complex tuyas, all of which have recognizable passage zones. Passage zones are mappable, stratigraphic boundaries that demarcate continuous transitions between subaqueous and subaerial volcanic deposits, for example, where subaerial lavas flow into a lake, river, or ocean (e.g., Jones and Nelson, 1970). Passage zones can also record highstands of intraglacial lakes formed during glaciovolcanic eruptions, and thus they can provide estimates  Note: Ages were calculated relative to 28.201 Ma for the Fish Canyon sanidine (Kuiper et al., 2008) using decay constants of Min et al. (2000). *Numbers in parentheses refer to volcanic center numbers on Figure 2.
of minimum ice thicknesses extant during the eruption (Jones, 1968;Russell et al., 2013Russell et al., , 2014cf. Smellie and Edwards, 2016). The other 10 newly analyzed deposits do not have obvious passage zones. At these locations, the entire volcanic pile apparently formed beneath the surface of a lake confined within an ice sheet. We interpret this to indicate that the ice surface must have been above the top of the present edifice. Although the 2σ errors for the new eruption ages prevent us from matching all 18 tuyas to single, unique marine isotope stages (MISs), the ages show that ice was present in northern British Columbia during 7 of the 13 glaciations since 1 Ma (MIS 4,6,8,16,18,22,24/26), and during at least five glaciations between 2.8 and 1 Ma.
Another five deposits are from subaerial eruptions. The oldest (4302 ± 24 ka) is a Pliocene lava on the surface of the Kawdy Plateau (Fig. 2). This age is consistent with the Northern Hemisphere record suggesting that regional ice sheets did not exist in the mid-Pliocene (Ehlers and Gibbard, 2008;Balco and Rovey, 2010). A 2226 ± 14 ka subaerial lava (MIS 85) that fills a large, U-shaped valley north of the Kawdy Plateau is the first record of relatively ice-free conditions after the onset of Cordilleran ice sheet formation recorded at 2.8 Ma as established here. A 1971 ± 15 ka subaerial lava shield is at a higher elevation, which prohibits even cirque glaciers at that time (MIS 75). A smaller lava plateau immediately northeast of Kawdy Plateau also lacks evidence for glaciovolcanism, documenting a period without regional ice sheets at 1657 ± 34 ka (MIS 57). The last clear evidence for the disappearance of regional ice before the Holocene is from subaerial deposits along Canyon Creek formed at 763 ± 13 ka (MIS 19). While we think it is unlikely that an ice sheet persisted across the northern Cordillera from 763 ka to 63 ka, we currently have no evidence strictly from the volcanic record for ice-absent conditions during that time period.
The geochronology results show that volcano-ice interactions have occurred in this volcanic field throughout the Pleistocene. A detailed discussion of the broader-scale implications of the tuya ages in a global context is beyond the scope of this work, but the new ages demonstrate how interspersed subaerial and glaciovolcanic deposits provide broad constraints on the timing and locations of ice. They also give constraints on posteruption modification times, discussed below.

PERIGLACIAL AND GLACIAL FEATURES OF THE TUYA-KAWDY REGION
The modern landscape of the Tuya-Kawdy region has abundant periglacial and glacial features, including well-sorted stone circles (Fig. 8A) and stone stripes (Fig. 8B). Given their size and position (<2 m in diameter/width), they may be evidence for past permafrost extent (Ballantyne, 2018) during MIS 1 and 2. Gabrielse (1998) noted these features as well as abandoned channels and isolated drumlins across the southern part of the Kawdy Plateau and into the Dease Lake 1:250,000 map sheet just south of the study area. Mathews (1947) reported an ice body on the northside of Kima 'Kho (his "Kawdy Mountain"), and two tuyas (Kima 'Kho and Horseshoe tuya) have unmodified moraines that may have formed during the Little Ice Age. Various other glacial features within the study area record the former presence of Pleistocene glaciers. While it seems likely that many, if not all, of these features formed toward the end of the Last Glacial Maximum (LGM; MIS 2), we know of no published work on exposure ages of these features in the study area.

Pleistocene Glacial Features
The Pleistocene glacial history in the Tuya-Kawdy area is recorded by small-scale striations on bedrock (Fig. 8C), megalineations tens of meters wide and hundreds of meters long (Fig.  8D), cirques (Fig. 8E), large U-shaped valleys (Fig. 8F), glacial erratics transported up to tens of kilometers (Fig. 8G), and isolated eskers (Fig. 8H) ( Table 2). Many of the localized features do not require the presence of an ice sheet (small-scale striations, cirques), but eskers in wide valleys, regionally transported erratics deposited on top of features more than 300 m above the surrounding terrain, and landscape-scale megagrooves and large U-shaped valleys are consistent with regional ice transport. Mathews (1947) reported that all of the tuyas he visited had glacial erratics.
Given the limited, regional-scale mapping in the study area, it is not surprising that deposits of glacial diamict have not been previously described. However, we identified a glacial diamict in the central part of the Kawdy Plateau, along the northern lower flank of Kima 'Kho tuya (Fig. 9). The deposit consists of poorly sorted and subrounded to subangular clasts (Figs. 9A-9C), some of them striated (Fig. 9C), and it is as thick as 25 m. The clasts include alkali olivine basalt, arkosic and quartz-rich sandstone, metagraywacke, and metarhyolite. While all of these rock types are present in the region, only the basalts could have been transported from local bedrock. The top of the deposit is friable material that is cemented at depth such that it forms low outcrops (Fig.  9D). Toward its southern end, the deposit has a planar contact with underlying tuff-breccia that is presumably part of Kima 'Kho tuya (1.8 Ma; Fig. 9D). We also used remote sensing to analyze the orientations of more than 800 mega-lineations on the central and southern Kawdy Plateau, to determine if they all record the same ice-flow directions (Fig. 10). Variations in these orientations could indicate either formation from different ice sheets or that ice-flow directions changed during the growth and decay of a single ice sheet. Where closely observed, the lineations are between 20 and 600 m long and generally less than 2 m high. Most lineations were found on the plateau surface (~1300 m a.s.l.), but they also mark the top and flanks of Horseshoe tuya (Fig. 10B, lowerright corner) and Tutsingale Mountain (Figs. 10C and 10D). The summit of Tutsingale Mountain is 420 m above the surrounding plateau surface, so the ice body that covered the mountain must have been significantly thicker. The summit lineations are oriented east-west and are obvious features preserved in the volcanic felsemeer that covers the summit (Fig. 9D). Statistical tests (one-way analysis of variance) confirmed that the lineations fall into three distinct subsets with orientations of 21°-55° (N = 23), 56°-86° (N = 751), and 87°-169° (N = 53; Fig. 10E). The dominant orientation was found across the entire Kawdy Plateau (Fig.  10E, i); the two subordinate orientations were only found on the flanks and summits of Tutsingale and Horseshoe tuyas (Fig. 10E, ii). It is possible that preservation of these lineations at higher elevations indicates that they formed from ice flow with a distinctly different direction than for LGM ice, or that ice-flow directions changed during the LGM.

Pleistocene Tuya Modification
In at least three locations, streamlined tuyas are also indicators of ice-flow directions and may mark the positions of paleoice streams (e.g., Eyles et al., 2018). The Tuya-Kawdy area presently has stream drainages that ultimately go to the Arctic Ocean, to the Bering Sea via the Yukon River, and to the Pacific Ocean ( Figs. 1 and 11). Given that modern rivers all flow through much larger, U-shaped valleys, we assume that basal ice flow during at least the most recent glacial stades was in these directions. In each of these drainage systems, at least one tuya has been shaped by posteruption ice flow. In the western part of the volcanic field, Klinkit Creek tuya sits in a small valley that drains into presentday Jennings River (Fig. 11B), which empties into the Bering Sea via the Teslin and Yukon Rivers. It has a central high summit dome (1490 m a.s.l.) at the upstream north end of the volcano, with a downstream-tailing basal platform (~1150 m a.s.l.). About 50 km west across the inferred drainage divide, Rancheria tuya has a summit plateau over 1700 m a.s.l. with an eastward tail at 1200 m a.s.l. (Figs. 11C and 11E). It sits within the watershed of Rancheria Creek, which flows into the Liard River via the Mc Kenzie River into the Arctic Ocean. The third example is South tuya, ~50 km south of Rancheria tuya and 8 km north of Tuya Lake (Figs. 11D and 11F). Its conical summit reaches to 1870 m a.s.l., overlooking a lower bench elongated to the south at ~1380 m a.s.l. It lies within the drainage of Tuya River, a major tributary of the Stikine River, which flows to the Pacific Ocean. The southwestern part of the field area has yet another regional outlet via the Nahlin and Taku drainages, but none of the tuyas in that part of the volcanic field is obviously streamlined, likely due to their position on top of the Kawdy Plateau.
The presence of cirques also records posteruption glacial modification (Table 4). Cirques are common throughout the field area (Fig. 2), with headwalls opening in all directions except south. The formation of cirques requires elevations high enough to maintain perennial snow for centuries at least. Given that the Pliocene-Pleistocene likely had many periods of warming and cooling coincident with marine isotope stages (Lisiecki and Raymo, 2005), it is possible that the Tuya-Kawdy area has experienced many periods of cirque glaciation. Yet, not all tuyas in the study area have cirques. This may be in part because many tuyas also lack lava caps that could facilitate cirque formation. Seven of the eight tuyas with obvious cirques are flat tuyas. We highlight four tuyas with different eruption ages that show distinct differences in their relative degrees of glacial modification (Fig. 12). Klinkit Lake tuya is essentially unmodified. It is a flat tuya formed ca. 0.12 Ma in the northern part of the field area, at the eastern end of a small mountain range just north of the large Jennings River valley (Figs. 2 and 12A). Its summit lava cap at 1520 m a.s.l. shows few signs of erosive modification or cirque development. Horseshoe tuya (moderately modified) has a single, well-developed cirque that appears to have cut almost halfway through its summit lava cap at 1790 m a.s.l. since it formed at ca. 2.07 Ma (Figs. 2 and 12B). It is located at a moderate elevation for the field area on the central part of the Kawdy Plateau that is flat except for the tuyas (see Fig. 3A). Mathews tuya (moderately modified) is located in the Cassiar Mountains just northeast of Tuya Lake between a large U-shaped valley to the south, and a smaller hanging valley to the north (Figs. 2 and 12C). While its summit elevation is similar to that of Horseshoe tuya (1795 m a.s.l.), it has almost lost its entire capping lava to a north-facing cirque, even though it is only about one third the age of Horseshoe tuya (Fig. 12C). Three-cirque tuya (heavily modified) has one of the highest summit elevations of all the tuyas studied (2050 m a.s.l.). Its heavily dissected lava cap hosts at least three cirques developed since 0.84 Ma (Figs. 2 and 12D). It also is within the Cassiar Mountains, on the western end of a ridge that ends in a large U-shaped valley. Most of the tuyas within the Cassiar Mountains have summit elevations below those of surrounding mountain peaks. Only one conical tuya in the study area shows any evidence of cirque development: Kawdy Mountain on the southwestern edge of the field area (see Fig. 10A). It also has a high elevation (>1900 m a.s.l.) and is old enough to have experienced several glaciations (0.96 Ma). However, it lacks any remnant capping lava. Since most of the other conical or linear tuyas do not show evidence of cirque development, perhaps Kawdy Mountain tuya originally had a lava cap that allowed cirque development but that has since been eroded away. Ash Mountain, the youngest (0.06 Ma; MIS 4) and highest (2100 m a.s.l.) tuya, has a conical morphology and records the highest minimum ice thickness in the study area, but it shows no sign of cirque development. It is also possible that either the lack of a subaerial lava cap or the lack of time impeded cirque development.

Ice-Thickness Estimates
While the LGM history of the Cordilleran ice sheet is well studied (e.g., Clague and Ward, 2011;Eyles et al., 2018), much less is known about the older history of the Cordilleran ice sheet (e.g., Barendregt and Irving, 1998;Froese et al., 2000). In the Tuya-Kawdy area, tuyas can constrain the presence, location, timing, and minimum ice thicknesses of the Cordilleran ice sheet for at least 17 out of 100 MISs over a span of 2.8 m.y. (Fig. 13)  showing major drainages within the field area, the locations of glaciovolcanic (red) and subaerial volcanic (green) deposits, and the locations of three ice-sculpted tuyas. Green shades on DEM are generally below 1100 m above sea level, and white shades are above 1500 m above sea level. (B) Klinkit Creek tuya is located on the east side of a U-shaped valley that drains to the south into the regional drainage of Jennings River. Ultimately, the water (and presumably former ice) flows into the Yukon River system, which empties into the Bering Sea. (C) Rancheria tuya is located on the north side of a large U-shaped valley that ultimately flows into the Liard River and the Arctic Ocean. (D) South tuya is on the south side of a small intramontane plateau, opposite Ash Mountain to the north. To the south, a large valley hosts Tuya Lake, which appears to be an ice-scoured lake basin. (E) View northward showing part of south flank of Rancheria tuya. The higher, west half of the edifice has a lava cap, but the east half appears to have been streamlined by glacier erosion. (F) View northward showing south side of South tuya. Deposits from this center end abruptly at its north edge but comprise a flank platform that extends 2 km south of the summit cone. from the eroded 2.8 Ma Meehaz Mountain on the western edge of the Kawdy Plateau (Figs. 2 and 13B) to little dissected 0.06 Ma Ash Mountain in the Cassiar Mountains (Fig. 13C). Massive subaerial lavas at Kima 'Kho directly overlie orange, subaqueous fragmental deposits in a spectacular passage zone (Fig.  13D). Generally, the subaqueous deposits contain pillows lavas as well as water-quenched, vitric ash-to block-sized grains that have been palagonitized, all directly overlain by subaerial lava (Fig. 13E). Passage zones can be used to reconstruct the elevations of intraglacial paleolakes: Where passage zones are present, they record minimum water levels of those lakes (Mathews, 1947;Jones and Nelson, 1970;Smellie, 2006;Edwards et al., 2009;Russell et al., 2013Russell et al., , 2014. The minimum paleolake surface also establishes a minimum height for the surrounding ice. The height estimate is considered minimum because the upper surface of many glaciers/ice sheets is permeable firn (50-100 m for modern glaciers; Cuffey and Paterson, 2010) that would allow water to drain (Mathews, 1947;Smellie, 2006). It is also likely that the level of the lake was much farther below the surface of the surrounding ice, but that is difficult to infer. At least one tuya, Kima 'Kho, has multiple passage zones that record changes in paleolake levels (Russell et al., 2013;Turnbull et al., 2016). The  Figure 12. Examples of the degree of glacial modification of individual tuyas. Each example shows an oblique view, a color-shaded digital elevation model, and a brief description. (A) Klinkit Lake tuya is a lava shield essentially unmodified by glaciation. (B) Horseshoe tuya has a summit lava shield into which one large cirque has been cut on its north side. (C) Mathews tuya has a summit lava shield that has a single cirque that has almost removed the entire lava cap. (D) Three-cirque tuya has been dissected by three cirques on the east, west, and north sides. Elevations are given in meters above sea level (a.s.l.).  level of the englacial lake would also have been controlled by the basal hydrology of the ice (lakes in frozen-based glaciers may not drain at all; those in polythermal or temperature glaciers may drain continuously) and by the local bedrock topography. Mathews (1947) envisioned "rocky spillways," the elevations of which controlled the levels and stability of the eruption-formed, intraglacial lakes. Tuyas without passage zones but consisting of lithofacies formed entirely beneath water also provide minimum icethickness estimates. At these volcanoes, the highest subaqueous deposit is the minimum estimate for the height of the enclosing lake (e.g., Ash Mountain; Fig. 13C). Entirely subaerial volcanic deposits are also critical paleoclimate proxies. Five subaerial lavas have been identified in the Tuya-Kawdy area. At least two, Jennings River and Canyon Creek, are flat-lying, valley-filling lavas that were emplaced into large U-shaped valleys (Figs. 13F and 13G). These occurrences record times when no regional, valley-filling ice sheet was present.
The compilation of tuya-constrained minimum ice thicknesses through time for the Tuya-Kawdy volcanic field shows that the Cordilleran ice sheet waxed and waned at least four times during the Pleistocene. While likely this happened much more frequently in concert with documented global glacial cycles, our study is one of the first to document relatively precise timing for Cordilleran ice sheet absence during the Pleistocene from the presence of low-elevation subaerial lava flows. The record also is consistent with an apparent Pliocene to late Pleistocene increase for minimum ice thicknesses in the Antarctic Peninsula (Smellie et al., 2008).

Ice-Extent Estimates
A secondary inference derived from minimum thickness estimates at tuyas is a minimum diameter for the enclosing ice. We think that each of the tuyas formed within a fully developed Cordillera ice sheet that spanned the entire field area and was continuous from the Coast Mountains to the Rocky Mountains. However, we do not yet have concrete evidence for a full-blown ice sheet independent of the volcanoes. Here, we estimated minimum ice extents as constrained directly by tuyas at specific locations and times to demonstrate that eruptions into local ice caps are not as feasible as eruptions into ice sheets. Nye (1952; see also Cuffey and Paterson, 2010) developed a paradigm for estimating the diameter of ice caps and ice sheets given ice-thickness estimates based on the assumption that those types of glaciers have a convex-up parabolic cross-sectional geometry. This approximation is useful for identifying places where a local ice cap (<50,000 km 2 in area) rather than regional ice sheet (>50,000 km 2 ) best explains the evidence for a paleo-ice body (e.g., Smellie and Edwards, 2016). For the Tuya-Kawdy tuyas, we used estimated minimum ice thicknesses (Fig. 13) and assumed that syneruption ice bodies were centered on each edifice. The estimated radii for the model ice bodies range from 2 to 28 km (Fig. 14A), with areas of ~12-2500 km 2 , consistent with local ice caps rather than a regional ice sheet (Fig. 14B). For example, the Chromite Creek glaciovolcanic deposit (Table 1; Fig. 2) is at an elevation of 1550 m a.s.l. So, when it formed (no age is available), the minimum local ice elevation must have been at least as low as 1550 m a.s.l. (Fig. 14C). If we assume regional ice covered all elevations above 1550 m a.s.l., ice would not have completely covered the landscape in northern British Columbia. The region potentially hosted numerous, disconnected small glaciers.
However, the physiographic positions of most of the ageconstrained tuyas are not generally consistent with "local" ice based on this same approach. The elevation of the Kawdy Plateau surface is ~1250 m a.s.l., which is at least 700 m below the elevation of mountain peaks and ridgelines in the Cassiar Mountains just to the north and the Coast Mountains farther west. If we use 1250 m a.s.l., those minimum ice elevations imply (e.g., Tuya Butte and seven others; Table 1; Fig. 14D) large ice caps and minimum ice elevations of 1050 m a.s.l. (most of the tuyas on the Kawdy Plateau require a regional ice sheet; Fig. 14E).

Pleistocene Landscape Evolution
Using this new glaciovolcanic record, we can speculate on the evolution of Pliocene-Pleistocene landscapes around the Tuya-Kawdy area of the North American Cordillera. With the modern topography and volcanic chronology as guides, we reconstructed the Tuya-Kawdy landscape in time (Fig. 15). Presumably, the preglacial topography in the western part of the area was still dominated by three elements: the relatively flat Kawdy Plateau to the south, several isolated massifs within the Cassiar Mountains, and large valleys separating the massifs (Fig. 15A). Mountain peaks were likely more rounded before the 2.8 m.y. of glacial sharpening. The first evidence for regionally extensive ice is recorded by Meehaz tuya at ca. 2.8 Ma; given its relatively low elevation on the Kawdy Plateau, possibly several of the mountain summits would have projected as nunataks above the regional ice surface (Fig.  15B). We note that existing models for formation of the southern part of the Cordilleran ice sheet show ice building up in the Coast Mountains in western British Columbia and in the Rocky Mountains in eastern British Columbia (Clague, 1989). Glaciers and ice caps originating within these highs would have partly flowed down into the interior of British Columbia to eventually coalesce into the Cordilleran ice sheet. While our new data and analysis document specific locations and times of ice presence, the minimum ice-thickness estimates from tuyas do not impose strict constraints on the slope of the ice-sheet surface in the Tuya-Kawdy area. The precision of the geochronology is not yet high enough to tie eruptions tightly to periods of ice-sheet waxing or waning. For our landscape reconstructions, we show the minimum ice-thickness constraints for a given time only to give a general sense of the minimum regional ice that might have been present. Meehaz tuya currently provides the oldest (2.8 Ma) direct constraint on ice-sheet formation in the Canadian Cordillera.
About 600 k.y. later, the subaerial Jennings River lava erupted into a wide valley that was at the time apparently ice free (Fig. 15C). Where exposed, this lava shows vertical, wide columnar joints indicative of subaerial emplacement. While it is possible that the surrounding mountains could have hosted small, local glaciers, the lack of evidence for lava-ice contact argues against a regional Cordillera ice sheet then. These two oldest ages are among the first direct links showing that at least one and possibly several ancestral Cordillera ice sheets formed and withered between ca. 2.8 and ca. 2.2 Ma. Less than 400 k.y. later (MIS 64/66), the complex tuya at Isspah Butte/Metah Mountain formed, possibly within an ice sheet thicker than the earlier one (400-500 m minimum). The regional ice implied by the passage zone at this tuya would likely have inundated the entire field area  Figure 14. Minimum constraints on syneruption radial ice distribution around each tuya based on lithostratigraphic estimates of minimum ice thickness (cf. Fig. 13) and the ice-cap model of Nye (1952). (A) Argon geochronology plotted against minimum estimated radii assuming eruptions were beneath an isolated ice cap. The model radii assume a parabolic decrease in ice thickness (e.g., Nye, 1952). Subaerial centers (open) plot at radii of zero. (B) Plan-view digital elevation model (DEM) of the field area showing the minimum allowable ice-cap areas. (C-E) Plan-view regional DEMs showing the surface cover of ice assuming complete ice cover above elevations of 1550 m above sea level (a.s.l.), 1250 m a.s.l., and 1050 m a.s.l. The approximate location of the regional map is shown in Figure 1.   Figure 15. Chronological snapshots of the minimum, syneruption ice distributions across the Tuya-Kawdy region implied by the minimum ice-thickness estimates derived from four separate volcanoes (red, glaciovolcanic; green, subaerial). Topographic profile is in the west half of the field area, from southwest (left) to northeast (see Fig. 2). Ice-sheet profile is shown with no surface elevation gradient because it is unconstrained at present. CIS-Cordilleran ice sheet; a.s.l.-above sea level. (A) Reconstructed topography before 2.8 Ma. The topography of the area is modified to approximate the preglacial profile. (B) Ice-sheet distribution at time of eruption of Meehaz Mountain tuya (2.8 Ma), implying a minimum number of nunataks in the Cassiar Mountains. (C) Interglacial period as recorded by the subaerial Jennings River lava formed at 2.2 Ma. (D) Glaciation recorded by the Isspah Butte-Metah Mountain complex tuya at 1.8 Ma. (E) Glaciation recorded by Klinkit Lake flat-topped tuya at 0.12 Ma, which has relatively low summit and effusive passage zone, and so records relatively low ice minimums. (F) Glaciation recorded during the Last Glacial Maximum assuming that erratics found on the top of Isspah Butte were deposited at the peak of marine isotope stage (MIS) 2. (G) Present-day topographic profile including volcanic centers erupted during the past 2.8 m.y. Peaks in the Cassiar Mountains are shown as "sharpened" by cirque erosion. (Fig. 15D). This possibly led to sharpening of any unglaciated peaks that protruded as nunataks above the ice.

Cassiar Mountains Kawdy Plateau
The last ice recorded by tuyas along this specific western transect was at 0.12 Ma, potentially during the transition from MIS 6 to MIS 5, when Klinkit Lake tuya erupted (Fig. 15E). Its apparent record of comparatively thin (see Table 4) regional ice could be recording Cordilleran ice sheet withering and disintegration. The modern landscape has not only been shaped by added volcanoes, but it also has been sharpened by repeated glaciations (Fig. 15F). Repeated formation and disintegration of the Cordilleran ice sheet have reduced deposit volumes of present-day tuyas. Presumably, most surviving periglacial features formed after the LGM, although differences in orientations of mega-lineations across the southern part of the Kawdy Plateau may signal multiple generations of ice-flow directions. If the Cordilleran ice sheet reached its maximum thickness during MIS 2, then the erratics found during this study and previous ones (e.g., Watson and Mathews, 1944;Gabrielse, 1998) on top of many tuyas show that they were all buried beneath ice at the LGM. While no single topographic profile through the field area adequately represents all glaciations recorded by tuyas described in this study, this region highlights the potential for using tuyas to inform the evolution of glacierized landscapes.

CONCLUSIONS
Abundant tuyas and glaciological markers in the Tuya-Kawdy region of northern British Columbia provide a clear record for the periodic existence of a northern Cordilleran ice sheet during the past ~2.8 m.y. Our 23 new 40 Ar/ 39 Ar ages show that regionally extensive ice was present at least sporadically from 2.8 to 0.06 Ma. The glaciovolcanoes broadly correlate with 14 of the globally recognized glacial marine isotope stages, and five subaerial deposits correlate with interglacial marine isotope stages. The 35 volcanic deposits of the Tuya Formation are dominated by glaciovolcanic lithofacies. Pillow lava and palagonitized tephra attest to water-dominated eruption environments in places where former ice is the only feasible water-confining agent. Tuya morphology and stratigraphic markers record minimum syneruption ice thicknesses. These can be used as important inputs for arguments to differentiate between local and regional ice bodies. The relative resistance of subaerial lava caps to glacial erosion reflects their physiographic positions and to a lesser degree their ages. The information from tuya ages, degrees of modification, and relative elevations can be combined to begin to assess Pleistocene landscape evolution beneath the many incarnations of the Cordilleran ice sheet.

ACKNOWLEDGMENTS
This work was initiated in 1995 by a grant from Lithoprobe to J.K. Russell, through which one of us (Edwards) was able to first visit the Tuya-Kawdy region. More recent fieldwork was largely funded by the Dickinson College Office of the Provost (2003 field season), the National Science Foundation (NSF) American Recovery and Reinvestment Act of 2009 (grant EAR-0910712 to Edwards, 2009Edwards, -2010, and the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery program (Russell). We thank the many students/collaborators who participated in fieldwork in this beautiful but remote region (C. Haynes, W. Kochtitzky, E. Lasher, C. Ryane, K. Simpson, K. Wetherell, M. Turnbull), as well as Jim and Sharon Reed for safe flying and logistical support. Reviews by M.T. Gudmundsson and J.C. Clague clarified the writing and organization of the manuscript, as did meticulous editing by R.B. Waitt.